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The  results  of  a head  injury  model  devel  ment  program  are  presented,  including  a des- 
cription of  the  resulting  model’s  features  and  its  capabilities  for  simulating  direct  and  indirect 
impact  forces.  The  model’s  validity  is  discussed  in  terms  of  level  of  confidence  and  verifica- 

Skull  bone  response  and  brain  response  are  presented  for  a variety  of  dynamic  simu 

lations.  Over  75  dynamic  and  static  computer  runs  have  been  executed  in  its  development. 
The  basic  features  of  the  model  are  described,  including  recognizable  skull  geometry,  linear 
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■elastic  and  linear  visco-elastic  behavior,  and  a capability  for  specifying  arbitrary  impact  loads 
and  boundary  conditions.  A special  modification  of  the  isoparametric  element  is  shown  to  be 
particularh  suited  to  simulation  of  the  dynamic  response  of  nearly  incompressible  brain 


A preprocessor  enables  automatic  mesh  generation  of  a skull  model  consistent  with  a 
prescribed  set  of  geometrical  data  supplied  by  the  user,  hither  complete  three-dimensional 
skulls  or  skulls  symmetrical  with  respect  to  the  midsagittal  plane  can  be  specified  in  the  mesh 
generation  process.*  Additionally , scale  factors  can  l>c  prescribed  which  modify  existing  skull 
meshes  and  achieve  parametric  control  on  si/c  and  shape.  A postprocessor  facilitates  the 
reduction  of  the  large  ambunt  of  data  that  is  typical  of  a head  impact  simulation.  The  scope 
anil  limitations  imposed  by  tbe  assumption  of  linearity  are  discussed.  The  results  demonstrate 
that  while  some  minor  changes  appear  indicated,  the  model  predictions  yield  useful  insight 
into  the  mechanical  causes  of  skull  and  brain  injury. 

Volume  1 of  this  report  also  contains  Appendix  A,  a clinical  description  of  head  injury. 
Volume  II  contains  Appendixes  B through  (i  covering  the  computer  programs  for  skull 
modeling. 
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The  results  of  a head  injury  model  development  program  arc  presented,  including  a des-  • 

cription  of  the  resulting  model’s  features  and  its  capabilities  for  simulating  direct  and  indirect  [ 

impact  forces.  Level  of  confidence  and  verification  of  the  model  is  discussed.  Skull  bone  and  brain  , 
responses  are  presented  for  a variety  of  dynamic  simulations.  Basic  features  of  the  model  are  des-  ' 
cribed,  including  recognizable  skull  geometry,  linear  elastic  and  linear  visco-elastic  behavior,  and  I 
capability  for  specifying  arbitrary  impact  loads  and  boundary  conditions.  A special  modification  of  • 
the  isoparametric  element  is  particularly  suited  to  simulation  of  the  dynamic  response  of  nearly  in  * 
compressible  brain  matter.  | 

A preprocessor  enables  automatic  mesh  generation  of  a skull  model  with  geometrical  data 
supplied  by  the  user.  Either  complete  three-dimensional  skulls  or  skulls  symmetrical  with  respect  ' 
to  the  midsagittal  plane  can  be  specified  for  mesh  generation.  Additionally,  scale  factors  can  modify^ 
existing  skull  meshes  and  achieve  parametric  control  on  si/e  and  shape.  A postprocessor  can  reduce  , 
the  large  amount  of  data  typical  of  a head  impact  simulation.  The  model  predictions  yield  useful  • 
insight  into  the  mechanical  causes  of  skull  and  brain  injury.  | 

Volume  I includes  a clinical  description  of  head  injury  (Appendix  A),  Volume  II  contains 
the  computer  programs  for  skull  modeling.  i 
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llj  weighting  factors  defined  in  Gaussian  quadrature  process 
k j j global  stiffness  matrix  coefficient  on  the  diagonal 
I t prescribed  impact  force  at  the  ir*'  degree  of  freedom 
I',-  displacement  at  ir*'  degree  of  freedom 
in  j lumped  mass  at  i1*'  degree  of  freedom 


xiii 


IN  I kODUC  I ION 


1 valuation  of  head  injury  received  from  direct  blows  is  necessary  for  development  of 
adequate  protection  systems  for  the  head.  The  need  for  biomechanical  models  that  predict 
head  injury  has  been  expressed  by  both  the  engineering  anil  the  medical  community:  e.g., 
to  aid  safetv  specialists  in  developing  protection  standards  for  vehicle  occupants;  to  aid 
engineers  in  choosing  militarv  aircraft  and  other  vehicle  designs  for  crew  and  occupant 
safetv  . to  improve  the  capability  of  neurosurgeons  in  managing  brain  trauma.  Irrespective 
ot  its  professional  application,  clinical-pathological  mechanical  correlation  of  head  injury  is 
what  is  sought.  Statistics  compiled  by  the  National  Safety  Council  reveal  that  the  vehicle 
accident  is  a major  cause  of  death  in  this  country.  Recent  reductions  in  maximum  speed 
limit  and  construction  of  safer  roads  in  this  country  are  improving  upon  the  highway  death 
rate  and  at  the  same  time  arc  making  the  crashworthiness  design  of  automobiles  more  fea- 
sible than  in  the  past,  f atalities  went  from  56,000  in  1973  to  46,000  in  1975'  . It  has  been 
estimated1  that  in  1975  the  Interstate  Highway  System  saved  4,400  lives.  Relative  reduc- 
tion in  fatalities  was  twice  as  great  on  those  roads  affected  by  the  speed  limit  reduction  to 
5 5 mph  as  on  those  roads  which  were  not  affected2 . The  goal  should  be  to  make  all  crashes 
survivablc.  With  regard  to  transportation  in  aircraft,  both  commercial  and  military,  the  goal 
cannot  lie  so  optimistic,  but  improvements  in  the  crashworthiness  design  of  aircraft  are  def- 
initely justifiable  and  warranted. 

Statistical  studies  of  injury  during  aircraft  crashes  have  shed  some  light  on  the  incidence 
ot  head  injury.  Two  difficult  questions  have  to  be  asked  of  the  data-.  “Were  the  crashes  in 
question  survivablc?”  and,  “In  fatal  instances,  could  death  be  attributed  to  head  injury?”  A 
survivablc  crash  has  been  defined  as  one  in  which  the  forces  transmitted  to  the  occupant 
through  his  seat,  restraint  system,  or  surrounding  components  do  not  exceed  human  toler- 
ance'’ . Though  if  is  the  best  we  have,  this  definition  is  only  preliminary  because  it  includes 
the  phrase  “human  tolerance.”  Certainly  one  of  the  goals  of  head  injury  research,  as  exem- 
plified by  the  model  studies  discussed  herein,  is  to  determine  the  levels  of  human  tolerance 
to  impact. 

Some  statistics  regarding  injuries  in  rotary  wing  aircraft  crashes  are  germane;  reference  3 
reports  on  2,546  helicopter  accidents  studied.  It  was  determined  at  the  outset  that  93%  of 
these  crashes  were  survivablc.  fhen  it  was  computed  that  40%  of  the  total  number  of  fatal- 
ities occurred  in  the  survivablc  group  (this  amounted  to  439  deaths  that  should  not  have 
occurred).  It  was  further  determined  that  2 3%  of  the  preventable  deaths  were  due  to  head 
injure  alone  and  that  29%  of  the  2,699  survivors  sustained  head  injuries. 


1 National  Safety  (Council  Accident  facts.  Chicago,  Illinois,  1976. 

2 

"National  Highway  ami  Iransportation  Safety  Agency.  NIII'SA  1)01  US  80171  5-.  I fleet  of  the  fuel  shortage  on  travel 
an»l  highway  safety,  by  I ( Cerrclli.  Washington,  I)  ( . Aug  1975. 

\ I*  Desjardins  Vehicle  crashworthiness,"  Numerical  and  Computer  Methods  in  Structural  Mechanics,  ed.  S.  J.  I enves, 
el  al  New  York.  Academic  Press,  1973,  pp  557  584 
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Impact  energy  associated  with  fixed  wing  aircraft  crashes  is  generally  expected  to  he 
higher  than  that  associated  with  rotary  wing  aircraft  crashes,  for  that  reason,  the  percent 
age  of  survivable  crashes  in  fixed  wing  aircraft  accidents  is  expected  to  he  less,  hut  crash- 
worthiness  of  these  aircraft  must  be  improved  as  well. 

I he  head  receives  basically  two  kinds  of  dynamic  loads  ( 1 > a direct  impact  involving 
hard,  localized,  contact-impact  forces  applied  directly  to  some  portion  of  the  skull  and  (2) 
an  indirect  impact  involving  forces  transmitted  from  the  neck  to  the  base  of  the  skull,  pos- 
sibly with  additional  nonlocali/ed,  soft  contacts  to  the  head.  I he  first  type  commonly  in- 
volves the  impact  of  heads  with  windshields,  instrument  panels,  sidewalls,  and  scat  hacks. 

The  second  occurs  when  the  occupant’s  restraint  system  (seat  belts,  air  bags,  etc.)  prevents 
his  reaching  surrounding  structures.  Space  for  occupants  is  always  at  a premium;  therefore, 
it  is  not  alwavs  feasible  to  remove  structural  components  from  the  occupants’  strike  zone4 
I he  only  alternative  is  to  design  the  occupants’  immediate  environment  so  that  if  contact 
occurs  the  effect  is  minimized  f.nergy-absorbing  structures  and  sometime  protective  helmets 
arc  essential  components  of  crashworthiness  design,  t heir  effectiveness  should  be  evaluated 
in  conjunction  with  a head  injury  model5,  along  with  other  surrogates  of  the  head. 

Approaches  to  Head  Injury  Prediction 

I he  prediction  of  head  injury  can  be  approached  in  many  ways-,  various  injury  or  s ver- 
ity indices  that  offer  attractive  expedient  solutions  have  been  proposed.  These  are  u.-ually 
based  on  data  derived  from  experiments  with  cadavers  and  human  volunteers.  This  data 
naturally  possess  limitations.  Often  little  agreement  exists  among  these  methods  of  injury 
piediction  because  each  method  chooscn  fits  the  experimental  data  with  a different  curve, 
but  the  primary  objection  to  these  methods  is  that  they,  by  definition,  must  either  simplify 
many  variable  relationships  or  fail.  As  a result  parameters  arc  often  contrived  that  have  no 
fundamental  basis  for  measurement,  but  which  are  critical  to  the  accuracy  of  the  method. 
Results,  therefore,  vary,  depending  upon  how  these  parameters  are  measured  or  estimated. 

In  any  approach  to  head  injury  prediction,  the  final  usefulness  of  the  derived  injury 
prediction  data  depends  on  how  conveniently  it  can  be  s\  nthesized  and  formatted  tor  assim- 
ilation. Primarily  because  of  the  computer,  this  convenience  requirement  docs  not  preclude 
a more  rigorous  approach  to  obtaining  useful  injury  prediction  information 
which  accounts  for  the  many  engineering  variable  relationships  present  in  the  problem, 
f urther,  if  the  more  rigorous  approach  is  successful,  more  is  obtained  than  an  estimate  ot 
injury  hazard;  it  also  provides  insight  into  the  clinical-pathological-mechanical  causes  of 
head  injury.  This  is  why  many  investigators  have  chosen  mathematical  modeling  as  an  ap- 
proach to  head  injury  prediction.  It  is  a more  instructive  approach  and  can  result  in  a useful 
and  convenient  data  base  of  head-injury  prediction. 


4 Army  Air  .Viol ■■  In y Research  ami  Development  I alioratory.  I K 71-22  Crash  survival  design  guide  I ort  I ustis,  Va, 

Ocl  1971 

’K  J Sac/alshi,  et  al  "A  critical  assessment  ol  the  use  of  non  human  responding  surrogates  tor  safety  system  evaluation," 
111  1‘roeectlmgs  of  chc  I wenticth  Srapp  . ar  ( rash  Conference,  SAP,  1976. 
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I his  research  was  to  develop  a computer  program,  based  on  the  laws  of  mechanics,  for 
simulating  structural  deformation  response  and  for  predicting  impact  injury  throughout  the 
head,  fhc  model  was  intended  for  use  in: 

(1)  understanding  injury  mechanisms 

(2)  defining  tolerance  envelopes 

(3)  predicting  dynamic  responses  to  specific  accidents 

14)  conducting  sensitivity  analyses 

1 he  specific  objective  was  to  compute  the  time  histories  of  displacement,  strain,  and  stress 
throughout  the  skull-brain  system  induced  by  arbitrary  head  impacts,  further,  a similar  pri- 
mate model  was  to  be  constructed  for  correlating  the  response  to  that  measured  in  animal 
impact  experiments. 

Scope 

Requirements  of  the  head-injury  model  (HIM)  included  three-dimensional  skull-brain 
geometry  simulation,  taking  full  advantage  of  symmetry  about  the  midsagittal  plane  in  the 
generation  of  the  geometrical  construction.  Neck-related  injury  was  not  included,  though 
arbitrary  boundary  conditions  at  the  base  of  the  skull  could  be  specified  to  approximate  the 
influence  of  the  neck  on  the  skull-brain  response  in  any  simulation.  Thus,  the  model  per- 
tained to  closed  brain  injury  and  skull  fracture  exclusive  of  any  mitigating  neck  influence. 
Arbitrary  initial  conditions  in  the  form  of  cither  initial  displacements  or  initial  velocities 
could  be  specified. 

A major  part  of  the  investigation  was  the  computation  of  impact  forces  existing  at  the 
interface  between  the  skull  model  and  arbitrary  simulated  targets.  I he  variation  in  contact 
torccs  and  areas  over  which  the  forces  act  (contact  area)  were  sought  as  a function  of  time. 
Non  linear  geometry  and  constitutive  behavior  were  to  be  included  in  the  contact  problem 
as  well  as  in  the  skull/brain  model. 

This  research  represents  3 years  of  study  under  contract  DOT-I IS-289-3-5501A  for  the 
Department  of  I ransportation,  Work  Request  N0020375WR00142  for  the  Naval  Aerospace 
Medical  Research  l aboratory  Detachment,  New  Orleans,  and  Work  Request  N0001 476WR- 
60083  tor  the  Office  of  Naval  Research. 
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I he  relatively  high  incidence  of  head  injuries  occurring  in  automobile  and  survivable 
aircraft  crashes,  combined  with  the  potential  seriousness  of  these  injuries,  has  been  respon- 
sible tor  the  preponderant  concentration  of  biomechanics  research  in  head  impact.  I .xperi- 
mental  research  has  sought  to  generate  data  in  support  of  head-injury  model  development, 
i ngineers  who  construct  models  must  be  cognizant  of  these  experimental  results.  Basically, 
models  must  be  capable  of  ( 1 ) accepting  as  input  the  mechanical  properties  of  the  head, 

(2)  including  both  experimentally  and  statistically  observed  mechanisms  of  injury,  anil  (3) 
simulating  live  primate  impact  tests  and  correlating  the  resulting  measured  data. 

Many  attempts  have  been  made  to  develop  biomechanical  models  for  predicting  mechan- 
ical phenomena  responsible  for  head  trauma.  Head  injury  studies  can  be  divided  into  two 
classes  primarily  experimental  and  primarily  theoretical.  I .xperimcntal  studies  utilize  an- 
thropomorphic dummies,  cadavers,  animals,  etc.,  and  require  expensive  and  time-consuming 
research  programs.  On  the  other  hand,  theoretical  studies  utilize  the  laws  of  mechanics  to 
predict  mechanical  responses  of  the  head  and  can  provide  an  efficient  tool  for  investigating 
head  injury  phenomena.  I undametallv,  only  two  kinds  of  mathematical  approaches  ana- 
lytical and  numerical  have  been  demonstrated  in  head  injury  modeling;  both  are  determin- 
istic. 

Analytical  Models 

The  analytical  models  appeared  first.  Kxamplcs  of  these  models  arc  one-dimensional 
fluid-tilled,  rigid  containers'’  and  two-dimensional,  axisymmctric,  fluid-filled,  elastic 
spheres  K I heir  solutions  were  usually  in  the  form  of  truncated  infinite  series,  and  assis- 
tance from  the  computer  was  required  to  obtain  data.  These  models  have  provided  useful 
insight  into  pressures  and  pressure  wave  propagations  within  a continuum  of  contained  com- 
pressible fluid.  However,  they  are  necessarily  limited  to  analytically  simple  modeling  para- 
meters, tor  example,  they  can  neither  simulate  recognizable  skull  geometry  nor  loading  his 
tories  and  boundary  conditions  for  head  impacts  typical  of  those  presented  in  vehicle 
accidents. 


I llayashi  Study  of  intracranial  pressure  caused  by  head  impact.”  Journal.  I acuity  of  I nginecring.  I’nivcrsity  of 
I okyo,  1 969, pp  30-59 

' \ I I ngin  ami  N King  l.iu.  Axisymmctric  response  of  a fluid-filled  spherical  shell  in  free  vibrations,*'  Journal  of 
Biomechanics,  vol.  3.  no.  1 . Jan  1970. 

HN  < l.ee  anti  S M Advani  transient  response  of  a sphere  to  torsional  loading  a head  injury  model  Math.  ' Bio 
sciences,  vol  r>.  1970,  pp  473  4Xf* 
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Numerical  Models 


I lie  belief  expressed  herein  is  that  only  the  numerical  models  have  the  potential  for  pro- 
viding all  three  of  the  necessary  and  basic  ingredients  stated  above  and  thus  offer  greater 
potential  for  reliable  prediction. 

f inite  Difference  Method.  The  finite  difference  technique  falls  within  this  category’; 
but,  in  spite  of  some  noteworthy  efforts,  this  technique  has  demonstrated  little  potential 
beyond  axisymmctric  simulation.  Reference  9 describes  a one-dimensional  finite  difference 
model  in  which  the  importance  of  including  the  various  layers  of  different  cranial  materials 
and  their  influence  on  dissipation  of  energy'  as  stress  waves  travers  the  layers  was  investi- 
gated. This  particular  work  w'as  then  extended  to  a two-dimensional,  axisymmetric 
sphere9 1 0 * ; others,  too,  have  discretized  axisymmctric  fluid-filled  shells  with  spherical1  1 and 
elliptical1 2 * 14 shapes.  Geometrical  zoning  in  the  finite  difference  technique  is  potentially  very 
flexible,  but  further  extensions  to  three-dimensional  head  injury  models  are  difficult  to  en- 
vision. l or  short  load  durations  typical  of  wave  propagation  regimes,  however,  the  finite 
difference  methods  are  attractive,  especially  when  nonhnear  behavior  and  axisymmetric 
skull  configurations  arc  considered. 

1 inite  Element  Method.  The  finite  element  method1 3 1 4 , on  the  other  hand,  is  an  es- 
tablished three-dimensional  modeling  technique  and  is,  therefore,  eminently  suited  for 
modeling  the  general  head  impact.  Therefore,  the  following  review  will  mention  only  those 
efforts  based  upon  the  finite  element  method.  A more  thorough  examination  is  made  in 
Reference  15.  Excellent  reviews  of  head  injury  models  based  upon  alternative  approaches 
can  be  found  in  Reference  16. 


9 

Nava!  Air  Development  ( enter.  Report  No  NAIX  -(  S 7113  Impact  analysis  of  the  skull-brain  system,  by  Stephen  I 
(•onion.  Warminster,  Pa,  Dec  1971. 

^Interim  Report  No.  NAIX  73065-40;  Analysis  of  head  impact,  by  Stephen  I..  (.onion.  Warminster.  I’a.  Apr  1973. 

1 Vl  V Benedict.  I II  Harris,  and  I)  U.  von  Rosenberg.  “An  analytical  investigation  of  the  cavitation  hypothesis  ot 
brain  damage. " Journal  of  Basic  I ngincering.  vol.  92,  Sep  1970.  pp  597-603. 

* “ Anthony  James  Oispino.  A dynamic  analysis  of  elastic  model  of  the  human  head.  M.  S thesis.  Department  ot  Mech- 
anical I ngincering,  University  of  Washington.  Seattle,  Wash,  1972 

* V)  ( /.ienkiewif;/.  The  finite  element  method  in  engineering  science  London,  I ngland,  Mc( iraw  Hill.  1971. 

1 4 

Richard  II  (.allaghcr  I- imte  element  analysis  I undamentals.  I nglewood  Clifts.  N.J ..  1975. 

Shugar,  I A..  “Simulating  and  Modeling  the  Human  Head's  Response  to  Impact,"  Aircraf  t Crashworthiness,  I ds 
K Sac/alski,  (•.  I . Singley  III,  W.  I).  Pilkcy.  and  R.  Huston,  University  Press  of  Virginia,  C’harlottsville.  197  5, 
pp  21  3-234. 

<} A l King  and  < < < hou  Mathematical  modeling,  simulation  and  experimental  testing  of  biomechanical  system 
crash  response."  paper  presented  at  I leventh  Annual  Meeting.  American  Institute  ot  Aeronautics  and  Astronautics.  Wash 
ington.  1)  C . I eb  1975.  ( AIAA  paper  no.  75-272) 
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The  first  suggestion  that  finite  element  techniques  could  be  advantageous  in  head  injury 
modeling  appeard  one  decade  ago1  7 , but  feasibility  of  a three-dimensional  model  at  that 
time  is  doubtful.  Axisymmetric  computer  codes  were  available,  however,  and  could  have 
easily  been  applied  toward  solving  the  fluid-filled,  spherical  shell  problem.  The  method  was 
not  applied  until  some  years  later1  N.  I his  application  was  significant  because  it  attempted 
to  deal  with  actual  skull  geometry,  and,  in  doing  so,  a primary  advantage  of  the  finite  ele- 
ment technique  in  head  injury  modeling  was  demonstrated  for  the  first  time.  The  analysis 
was  limited  to  an  empty  skull  and  static  loading  but  did  conclude  that  the  sandwich  nature 
of  cranial  bone  should  be  modeled  as  a layered  structure. 

Axisymmetric  f inite  l.lement  Model.  An  axisymmetric  finite  element  model  was  con- 
structed to  demonstrate  the  effectiveness  of  the  method  as  a potential  head  injury  investi- 
gatory technique1  v . Spherical  configurations  were  assumed  to  be  useful  head  injury 
models;  and,  if  the  finite  element  model  correlated  well  with  experimental  data  obtained 
from  an  instrumented  spherical  aluminum  shell,  the  model  was  assumed  to  be  reliable. 
Computed  strains  were  about  20%  higher  than  measured  strains  on  the  aluminum  surface. 

A subsequent  attempt  at  axisymmetric,  finite  element  modeling  of  the  skull  was  con- 
ducted with  hi  ore  success20.  The  axisymmetric  configuration  was  more  justifiable  when  the 
analysis  of  protective  helmets  was  undertaken.  Rotational  symmetry  is  a convenient  and 
reasonably  accurate  configuration  for  evaluating  a variety  of  helmet  parameters,  even 
though  the  head  itself  is  not  rotationally  symmetric,  furthermore,  when  a two-dimensional 
analysis  can  accomplish  the  objective  sufficiently,  it  is  to  be  preferred  over  a fully  three- 
dimensional  analysis,  f indings  regarding  helmets  will  not  be  discussed  in  this  report,  but 
such  work  does  offer  some  significant  points  regarding  head  injury  modeling  with  the  finite 
element  method.  An  inviscid  fluid  was  reported  as  being  employed  for  brain  material  char- 
acterization during  analyses  with  one  axisymmetric  code.  A second  computer  code  capable 
of  specifying  linear  viscoelastic  properties,  but  which  apparently  could  not  also  specify  the 
fluidity  alluded  to  earlier,  was  en. ployed  in  the  simulation  of  an  empty  spherical  shell.  The 
consitutive  relation  for  the  shell  was  an  exponentially  decaying  relaxation  modulus  devel- 
oped experimentally  from  impact  load  data  on  bone  samples.  The  experimental  data  from 
which  the  relaxation  modulus  was  derived  were  characterized  by  load  durations  in  the 
microsecond  range  and  were  applicable  to  impact  loads  with  durations  of  approximately 
50  /a sec  and  less.  Rather  good  correlation  was  shown  in  a comparison  of  computed  and 
measured  data  for  the  axisymmetric  configuration.  The  axisymmetric  trest  setup  displayed 
support  conditions  which  were  rotationally  symmetric.  The  spherical  aluminum  shells  were 
covered  with  a nested  set  of  spherical  shell  caps  which  simulated  the  helmets. 


* 7W  Goldsmith  I he  physical  processes  producing  head  injury,”  in  Proceedings  of  the  I lead  Injury  Conference,  I ippin 
cott,  Philadelphia,  Pa.  1966,  pp  350-382. 
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Office  of  Naval  Research.  Technical  Report  No.  8 Plastic  analysis  of  a skull,  by  ( II.  Hardy  and  P.  V.  Marcal.  Wash 
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V II  Kenner  and  W Goldsmith.  “Dynamic  loading  of  a fluid  filled  spherical  shell.”  International  Journal  of  Mechanical 
Silences,  vol  14.  1972.  pp  557  568. 

~M  | B Khali,  W Goldsmith,  and  J I Sackman.  “Impact  on  a model  head-helmet  system,”  International  Journal  of 
Mechanical  Sciences,  vol.  16,  1974,  pp  609-625. 
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It  is  significant  to  note  that  when  the  axisymmctric  finite  element  results  were  com- 
pared with  experimental  data  obtained  from  a cadaver  head  fitted  with  a helmet,  the  corre- 
lation deteriorated . I he  reasons  stated  were  that  the  boundary  condition  or  geometry  of 
the  cadaver  experiments  were  not  faithfully  duplicated  and  that  correct  material  behavior 
was  not  sufficiently  simulated. 

Spherical  and  ellipsoidal  models  were  constructed  in  a more  recent  two-dimensional 
finite  clement  analysis  study2 1 to  devise  a better  model  one  that  would  reveal  the  etiology 
of  head  injurs  due  to  impact.  I he  ellipsoidal  model  seemed  to  be  a slight  improvement  over 
the  spherical  model;  nevertheless,  the  results  were  stated  to  be  inconclusive.  I he  computed 
intracranial  response  was  based  on  a viscoelastic  core  material,  which  was  specified  in  the 
form  of  a Voight  solid  for  the  shear  moduli;  the  bulk  moduli  were  specified  as  elastic  con- 
stants. An  uncertainty  was  expressed  concerning  the  viscous  component  constant. 

Another  interesting  feature  of  this  finite  element  study  was  the  manner  in  which  the 
temporal  integration  of  the  system  equations  of  motion  was  carried  out.  I he  iterative 
method  is  not  commonly  employed  in  the  mainstream  of  dynamic  finite  clement  analysis. 

It  holds  some  promise  for  large  nonlinear  systems,  but  currently  suffers  from  two  objec- 
tions- - f irst,  convergence  to  required  accuracy  has  been  found  to  be  slow  for  general 
applications,  making  them  inferior  to  direct  integration  methods.  Second,  it  is  impossible 
to  repeat  solutions  lor  a dillercnt  set  ol  loading  conditions  without  once  again  sustaining 
the  cost  of  the  iterative  integration  process.  The  result  is  that  the  iterative  method  can  be 
more  expensive  than  the  more  common  direct  integration  methods  for  linear  systems.  I he 
study  does,  once  again,  demonstrate  the  utility  of  the  finite  clement  method  in  head  injury; 
m this  case  the  effects  of  geometry  and  duration  of  impact  were  easily  examined.  Also  of 
importance  is  the  choice  of  a computer  code  that  allows  for  viscoelastic  behavior. 

A three-dimensional  finite  element  model  of  the  brain  applicable  to  prescribed  transla- 
tional acceleration  of  the  skull  symmetrical  about  the  midsagittal  plane  was  constructed2  3 . 

A rigid  skull  was  assumed  and,  therefore,  was  applicable  to  brain  injury  in  those  cases  where 
skull  deformation  does  not  contribute  to  the  injury.  Obviously  the  model  cannot  directly 
relate  to  skull  fracture.  More  recently  the  model  has  been  extended  to  account  for  rota- 
tion2 4 This  work  demonstrates  the  detail  which  can  be  achieved  in  simulating  geometry. 


II  S < han  "Mathematical  model  lor  dosed  head  impact,"  in  Proceedings  of  the  eighteenth  Stapp  Car  Crash  Confer- 
ence. Ann  Arlmr.  Mich.  Dec  1974,  pp  557-578. 

""I  rnst  Schrcm  Computer  implementation  of  the  finite  element  procedure,"  Numerical  and  Computer  Methods  in 
Structural  Mechanics,  ed.  S J.  lenves,  et  al  New  York,  Academic  Press,  1975.  pp  79  1 17. 

( 1 A dynamic  finite  element  model  of  the  human  brain.  Ph.  D.  thesis,  Department  of  engineering,  University 

Of  California  at  Cos  Angeles.  Cos  Angeles,  Calif,  1974. 

( C Ward  and  R B Thompson.  “The  development  of  a detailed  finite  element  brain  model,"  in  Proceedings  of  the 

Nineteenth  Stapp  Car  Crash  Conference,  San  Diego,  Calif,  Nov  1975.  pp  M1-A7M 
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Various  intracranial  structures  were  discretized  in  an  extremely  painstaking  fashion,  lore- 
most  among  them  were  the  falx  and  tentorium  membranes  that  divide  the  cranial  vault  into 
separate  cavities.  In  this  work  a Poisson’s  ratio  of  0.48  and  a Young’s  modulus  of  10  psi 
were  specified  for  brain  material.  At  first  glance,  these  numbers  appear  insignificantly  dif- 
ferent from  published  values  but  in  fact,  result  in  a bulk  modulus  of  approximately  80  psi . 

I bis  value  represented  a compressible  brain  material  characterization  whose  resistance  to 
volumetric  change  was  over  three  orders  of  magnitude  less  than  the  reported  value  of 
305,000  psi’ ' . No  provisions  were  made  for  the  inclusion  of  viscoelastic  material  charac- 
terization. A modal  analysis  or  eigenvalue  approach  was  employed  to  solve  the  equations  of 
motion. 


*"  "'West  Virginia  University,  Biomechanics  Laboratories,  Department  of  I ami  AM.  I lead  ui|iiry  model  construction  pro 
gram  data  compilation  and  review  Morgantown,  W.  Va.,  Jun  l‘>71  (Contract  No.  IMI-43-67  1 137) 
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The  problem  of  developing  a head  injury  model  is  of  a complexity  different  from  that 
usually  encountered  by  engineers.  A cautious  and  deliberate  approach  was  therefore  taken 
Several  research  groups  were  ultimately  involved  in  the  development  of  the  head  injur) 
model.  I o facilitate  a discussion  of  their  participation  a schematic  diagram  of  the  head 
injure  model  development  program  is  presented  in  Figure  2-1 . The  sponsoring  agencies 
\l  I I SA,  N AMR  I . and  ONR*  gave  CHI.  overall  responsibility  for  carrying  out  the  various 
research  tasks,  each  of  which  is  designated  with  a box  and  located  within  the  dashed  line. 

I hese  tasks  were  either  performed  in-house,  or  a subcontract  was  let  to  an  outside  research 
group. 

Simulation  of  the  cranial  anatomy  as  closely  as  possible  was  first  necessary.  I his  re- 
quirement superseded  all  other  requirements  in  the  model  development  because  it  distin- 
guished this  head  injury  model  effort  from  all  previous  efforts.  However,  reasonable  limits 
exist  for  geometrical  simulation  in  any  model  development.  A subcontract  was  thus  ini- 
tiated to  delineate  those  areas  of  the  skull/brain  anatomy  clinically  known  to  have  a higher 
frequency  of  involvement  with  head  trauma.  More  attention  was  to  be  given  those  areas 
during  the  process  of  discretization.  Clinical  descriptions  of  head  injury  mechanisms  were 
obtained  so  that,  where  possible,  the  model  could  be  constructed  to  accommodate  such 
mechanisms.  The  Los  Angeles  County/University  of  Southern  California  (LAC/USC)  Medi- 
cal Center  furnished  this  information ; their  report  is  included  as  Appendix  A. 

Developing  a three-dimensional  discretization  proved  to  be  a difficult  task  and  required 
more  effort  than  originally  expected.  Manual  discretizing  of  the  skull  continuum  was  ini- 
tially attempted  in  an  effort  to  be  extremely  responsive  to  geometrical  detail.  1 lowever, 
measuring  difficulties  and  uncontrollable  element  aspect  ratios  prevented  success.  Instead, 
automated  procedures  proved  to  be  the  optimal  method,  although  some  accuracy  may  have 
been  compromised  in  the  simulation  of  geometry. 

The  skull/brain  module  was  to  be  actually  the  end  product  of  the  discretization  process 
and  entails  more  than  the  automatic  mesh  generation  scheme.  For  example,  it  includes 
parametric  control  on  size  and  shape  of  the  skull  being  discretized,  a mesh  checking  scheme, 
mass  distribution  computation  capability  (skull  inertia  tensor),  and  several  default  options, 
all  of  which  attempt  to  facilitate  the  user  in  his  effort  to  construct  a three-dimensional 
discretization. 

I he  deliberateness  of  the  approach  was  perhaps  nowhere  else  more  apparent  than  in  the 
construction  of  successive  models.  I'hc  development  of  one-dimensional  (actually  a pris- 
matic assemblage  of  three-dimensional  elements)  models  and  the  gradual  increase  in  the 
complexity  of  two-dimensional,  axisymmetric  and  plane  strain  models  instead  of  concen- 
trating on  a three-dimensional  model  from  the  start,  was  the  intent  ol  the  successive  model 


•National  Highway  I raffu  Safety  Administration,  Naval  Aerospace  Medical  Research  I aboratory  Detachment,  Michoud 
Station,  and  the  Office  of  Naval  Research. 
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approach.  In  spite  of  their  obvious  spatial  differences  the  models  were  closely  related  nu- 
merically, and  because  of  this  the  successive  model  approach  is  especially  effective  in  an 
orderly  development  of  the  final  three-dimensional  model.  A relatively  new,  general  pur- 
pose, finite  element  analysis  program  (IEAP)26 , developed  at  the  University  of  California 
at  Berkeley  (UCB),  was  employed  in  the  successive  model  approach.  In  the  process,  each 
model  was  composed  of  elements  from  the  same  family  of  isoparametric  finite  element  for- 
mulations, and  each  was  integrated  in  the  time  domain  by  the  same  method;  i.e.,  the  same 
implicit  integration  operator  was  employed  throughout  the  development.  Therefore,  infor- 
mation and  experience  gained  from  exercising  the  preliminary  models  were  ultimately  appli- 
cable to  the  numerical  behavior  of  the  three-dimensional  model. 

It  w as  recognized  early  that  the  basic  I-  KAP  code  would  have  to  be  modified  in  two 
areas  to  facilitate  development  of  the  three-dimensional  head  injury  model.  A linear  visco- 
elastic capability  was  thought  to  be  highly  desirable  in  view  of  the  large  volume  of  literature 
on  viscoelastic  properties  of  biomaterial27  3 2.  A subcontract  was  let  to  the  University  of 
California  to  incorporate  into  the  code  the  capability  of  prescribing  linear  viscoelastic  mech- 
anical properties  for  the  skull/brain  materials. 

It  was  also  recognized  that  a finite  element  head  injury  model  would  be  relatively  ex- 
pensive to  operate  parametrically.  Computer  costs  associated  with  dynamic,  three-dimen- 
sional models  of  any  kind  are  known  to  be  high;  but,  add  to  that,  construction  of  a model 
where  geometrical  detail  is  of  high  priority  and  costs  arc  pushed  even  higher.  To  effect 
modifications  which  would  improve  the  I-  l-.AP  code’s  efficiency,  a subcontract  was  let  to 
I ngineering/ Analysis  (E/A)  Corporation  of  Kedondo  Beach,  California.  A summary  of  the 
modifications  is  presented  in  Appendix  B (Volume  II).  K/A  also  provided  a bandwidth 
minimi/er  program  which  optimally  reorders  the  nodal  point  numbering  system  for  any 
skull/brain  discretization  prior  to  its  use  as  input  data  for  the  head  injury  model  code. 


J’l  niurMiy  of  California,  Department  of  Civil  I ngineering.  I mite  element  analysis  program,  l>y  l<  I Taylor  Berkeley, 
< alit 


27 

J VV  Pugh,  ct  al  I iastu.  an. I viscoelastic  properties  of  trabecular  bone  impendence  on  structure."  journal  of  Bio 
mechanics.  vol  6.  no.  5,  sep  1973 
2H 

N < lung  Stress-strain  history  relations  of  soft  tissues  in  simple  elongation,  " Biomechanics,  ed  \ ( Tung,  \ 
Pcrrone,  and  M Anlikcr  Inglewood  ( liffs,  NJ.  Prentice-llall,  1972  pp  IH1-20K 
29 

K t lennyson.K  I wert.  and  V Niranjan  Dynamic  viscoelastic  response  of  bone,"  I xperimental  Mechanics, 

Nov  19  72 

U||  ( Wang  and  A S Wineman  "A  mathematical  model  tor  the  determination  of  viscoelastic  behavior  of  brain  in  vivo 
I Oscillatory  response,"  Journal  of  Biomechanics,  vol  5.  1972.  pp  431  446. 


iv  72 

J2m 

1 970 


I ruong  Visco-elastic  propagation  ol  longitudinal  waves  in  skeletal  muscle,"  Journal  of  Biomechanics,  vol  s. 
pp  I lO 

Cal  ford  and  J II  Me  I Ihaney.  A viscoelastic  study  ol  scalp,  brain,  and  dura,"  Journal  of  Biomechanics,  vol  3. 

pp  21  I 221 


A comprehensive  program  ot  head  injury  model  validation  was  developed  in  consonance 
with  Nl  1 1 SA.  ( I 1 , though  not  involved  directly  with  the  experimental  work  itself,  helped 
establish  the  validation  test  plan  requirements  for  experimental  data  so  that  the  kinds  of  ex- 
perimental tests  conducted  could  also  be  simulated.  Nl!  I SA  contracted  with  the  Highway 
Safety  Research  Institute  (HSRI),  University  of  Michigan,  to  design  and  conduct  experi- 
mental tests,  the  data  from  which  would  constitute  the  primary  means  for  verifying  the  ac- 
curacy of  the  head  injury  model.  The  types  of  experimental  tests  included  ( 1 ) static  load 
deflection  measurements  of  dry  rhesus  monkey  skulls,  (2)  dynamic  tests  w ith  strain  gages 
affixed  to  the  same  skulls,  and  (3)  impact  tests  on  anesthesi/.ed  rhesus  monkeys  instru- 
mented with  both  strain  gages  and  pressure  transducers. 

After  completion  of  the  validation  effort  a sensitivity  analysis  or  parameter  study  was 
conducted  wherein  a series  of  dynamic  simulations  were  made  using  the  displacement 
boundary  conditions  at  the  base  of  the  skull  as  the  parameters.  Among  the  many  alternative 
parameter  studies  that  could  have  been  conducted,  this  particular  study  was  found  to  be 
more  useful  because  the  head  injury  model  does  not  possess  a neck  discretization  and,  there- 
fore, cannot  simulate  directly  the  influence  of  the  neck  on  the  base  of  the  skull  during  im- 
pact. However,  by  investigating  the  sensitivity  of  the  skull  bone  structural  response  and  the 
intracranial  pressure  response  to  changes  in  the  prescribed  displacement  boundary  condition 
at  the  skull  base,  the  limits  of  the  neck’s  influence  can  be  established. 

I'o  satisfy  the  requirement  for  a nonlinear  head  injury  model  the  linear  IT.AP  code  had 
to  be  extended  to  the  nonlinear  regime.  I hough  some  provision  was  made  for  nonlinear 
behavior  in  the  original  code  it  was  decided  to  proceed  with  a linear  model  development  and 
to  subcontract  the  nonlinear  code  development.  This  work  was  eventually  accomplished  by 
the  original  I TAP  code  author  at  UCB. 

I wo  additional  requirements  of  the  head  injury  model  were  associated  with  the  nonlin- 
ear development.  A nearly  incompressible,  nonlinear  finite  element  had  to  be  added  to  the 
element  library  to  model  the  intracranial  contents.  Also,  a general  contact/impact  force 
computation  capability  needed  to  be  developed  to  predict  accurately  the  stress  histories 
existing  at  the  interface  between  the  skull  and  an  arbitrary  target  Individually  or  collec- 
tively, these  requirements  suggested  a basic  research  project,  as  contrasted  with  the  linear 
model  development,  and  as  such  were  ultimately  handled  separately  by  subcontract  to  UCB. 
Much  progress  was  reported  in  a series  of  documents  ’7  , and  the  final  code  implementa- 
tion of  the  developed  theoretical  concept  has  been  completed  in  two  dimensions. 

< i vi  I I ngm<  ermg  I aboratorv  < K75.<K»7  I*  in  it  c*  element  formulation  anti  solution  ol  contract -impact  problems  in  con 
ionium  mechanics,  by  I I K Hughes,  K l I ay  I or  anti  j I Sackman  Cmversny  ot  California,  Berkeley,  Calit.  Way 
1974  (SI  SM  Report  No.  74  K) 

R75  (M1H  I mite  element  formulation  and  solution  ot  contact  impact  problems  in  continuum  mechanics  Part  II,  by 
I I K Hughes,  |<  | I avlor,  and  | I Sackman.  I niversit \ ol  ( aliformu.  Berkeley.  ( aid,  Jan  1975  (SI  SM  Report  No 
75-3 ) 

( K77.001  I inite  clement  formulation  and  solution  of  contact -impact  problems  in  continuum  mechanics  Part  III,  by 
I | R Hughes.  R | I aylor,  and  I I Sackman.  I 'niversit  y of  (.ahtonua,  Berkeley,  < aid.  Jul  1975  < S I SM  Report  No. 

75  7) 

( ivil  I nginct  ring  I ahorator)  < K77.0U2  I inite  element  formulation  and  solution  ot  contact  impact  problems  m con 
tmuuni  mechanics,  Part  IV  by  I I R Hughes,  R | l avlor,  J R Sackman.  and  W Kanoknukulchai  l niversity  ot  ( al- 
iforma  Berkeley . Calif,  1976 

l I R Hughes,  ei  al  I inne  element  formulanon  and  solution  ot  a class  ot  contact  impact  problems  in  continuum 
mechanics,"  paper  presented  .11  I lord  Conference,  Structural  Mechanics  Reactor  technology.  I ondon,  1975. 


3.  I III  OKI  I 1<  \l  ( ONSIDI  K \ I IONS 


In  tins  section  the  gov  erning  matrix  equations  of  motion  lor  dynamic  linear  viscoclas- 
tieity  are  developed  beginning  with  the  theorem  of  virtual  work.  I hen,  before  their  solution 
is  described,  two  separate  discussions  are  given  which  describe:  first,  how  the  HIM  code- 
deals  efficiently  with  viscoelastic  history  vectors,  and  second,  how  the  nearly  incompressible 
property  of  brain  material  is  handled. 

Principle  of  Virtual  Work  and  I quations  of  Motion 

I he  principle  of  virtual  work  states: 

If  a structure  is  in  ‘‘equilibrium"  under  a set  of  external  forces  and  if  the 
structure  is  given  a virtual  displacement  consistent  with  the  constraints  of 
the  structure,  then  the  external  work  done  is  equal  to  the  internal  virtual 
work  done. 

I he  word  equilibrium  is  placed  in  quotes  because  it  can  refer  to  either  static  or  dynamic 
equilibrium.  Strictly  speaking,  with  the  latter  form  we  arc  deal ing  with  not  the  original 
principle  of  virtual  work,  but  with  D’Alemberts’  extension  of  the  original  version  which  in- 
cludes the  inertial  forces.  I he  principle  of  v irtual  work  applies  to  viscoelastic  constitutive 
relations  as  well  as  linear  elastic  constitutive  relations.  In  this  presentation  the  inertial  terms 
and  the  constitutive  relations  are  included  in  generalized  form  for  completeness  because  the 
head  injury  model  possesses  dynamic  capabilities  and  both  linear  elastic  and  linear  visco- 
elastic materials  capability. 

The  virtual  work  principle  may  be  expressed  for  a body  as 

8W|  6VVj.  0 (3-1) 

where  6W|  and  6Wj.  are  the  internal  and  external  virtual  work  expressions,  respectively. 

They  are  defined  as  follows 


6W|  = ,|y/>{(  } 1 f„}dV 


(3-2) 


5W,  - /vft{d}  1 {f}dV  + /sft{d}  1 {p}dS 


(3  3 ) 


where  {tlandio}  arc  the  strain  and  stress  vectors,  respectively 
V is  the  volume  of  the  body 
{d}thc  displacement  vector 

{l  [is  the  body  force  vector  containing  the  inertia  forces 
i p}  the  external  surface  traction  vector  prescribed  over  the  surface  S 
6(  ) the  f irst  variation  of  ( ). 


At  this  point  it  is  convenient  to  discretize  liquation  3 1 in  developing  the  finite  element 
formulation.  However,  this  can  be  done  in  either  of  two  ways.  1 he  first  is  to  discretize 
I quation  3-1  exactly  as  it  is;  i.c.,  pertaining  to  the  whole  body.  An  alternative  way  is  to  re- 
write liquation  3-1  for  a single  finite  element  anti  presume,  for  the  time  being,  that  the 
body’s  virtual  work  can  be  expressed  legitimately  as  a summation  of  all  elements  of  the  vir- 
tual work  contribution  from  each  element.  I he  latter  approach  is  convenient  in  this  case 
because  it  takes  full  advantage  of  the  finite  element  method’s  ability  to  account  for  the  pre- 
scription of  different  constitutive  relations  within  the  skull/brain  continuum. 

I'he  basic  finite  element  chosen  in  constructing  the  head  injury  model  is  the  eight-node, 
isoparametric  element  and  is  shown  in  figure  3-1  together  with  its  assumed  linear  displace- 
ment function. 

The  X-Y-Z  axes  form  an  inertially  fixed  reference  system.  I'he  displacements  u,  v,  and 
w are  defined  with  respect  to  this  system.  The  ij-rpf  axes  define  the  local  element  coordi- 
nate system.  The  displacements  within  an  element  and  on  its  boundaries  are  functions  of 
the  local  coordinates  in  the  form  of  shape  functions  N:  and  the  nodal  point  displacements 
at  the  corners  Uj,  vq,  and  Wj. 

A general  definition  of  strain  is  given  by  Green’s  strain  tensor  for  which  tx  and  exv  are 
two  examples  of  the  six  independent  components  of{t} . fhese  two  components  are 
defined  as 


(3-4) 


(3-5) 


I or  a linear  analysis,  it  is  assumed  that  the  strains  are  small  and  their  products  (second-order 
terms)  may  be  neglected.  I o this  small  strain  assumption  is  added  a small  displacement  as- 
sumption Practically,  this  means  that  the  geometry  of  the  elements  remain  basically  un- 
changed during  the  loading  process  and  that  first  order,  inf  initesimal,  linear  strain  approxi- 
mations can  be  used.  In  this  way  Equations  3-3  and  3-5  reduce  to 


Ax 


. 3 u 


9 x 


(3-6) 
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and 


. 1 /i)u  3v 
lx>'  2 yi)y  + 3x 

Because  the  displacement  components  are  functions  of  the  local  coordinates,  the  partial 
differentiation  indicated  must  be  carried  out  using  the  chain  rule.  I he  differentiation  will 
result  in  a strain-displacement  matrix  [B]  whose  components  are  functions  of  the  local  co- 
ordinates; i.e.. 


H = IK<S.r?,f)l 

where  the  generalized  nodal  point  displacements  are 


(3-8) 
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(3  9) 


We  may  now  express  the  body  force  vector  for  the  element  in  terms  of  the  corresponding 
generalized  acceleration  as 


|f|  = (Pi  Jdf  Ip]  IN]  |_5_ 
anil  the  external  traction  orce  vector  as 


(.10) 


i’xL 

JV, 

Pzi 


(3-11) 


I quation  3-1 , the  virtual  work  equation,  may  now  be  written  using  the  generalized  dis 
placement  and  the  strain  displacement  relationship  for  the  eight-node  element  as  follows 


fn  I»»r''j«|dV  + /„ 

V 


6iij 

.^i  | IN]  1 (p]  (N|  ( vi 

(>W;  \ I W, 


dV  fp 
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while  the  superscript  n re  ter  s to  the  volume  or  surface  area  ot  the  n^'  element  o!  the  body. 

At  this  point  w e can  either  specify  a line  tr  elastic  or  linear  viscoelastic  r<  itionship  be 
tween  stress  Jo}  and  strain  j < | 

1 he  linear  elastic  constitutive  relationshi  > between  stress  and  strain  is  g n by 


It  Kquation  3-1  3 is  substituted  into  liquation  3-12  and  the  arbitrary  virtual  displacements 
factored  out,  we  have 
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I he  stillness  matrix  |k|  ” ami  mass  matrix  [m] 11  of  the  n^1  element  are  defined  as 

Ik]"  = f n [B]  1 [1)1  I B]  dV 
V 

and 


[m  1 11  = fn  [N  | 1 [p  1 [N ) dV 
V 

Carrying  out  the  integration  in  local  coordinate  space  vve  first  make  use  of  the  relationship 
between  elemental  volumes  in  the  two  coordinate  svstems 


dV  = dX  dY  dZ  | J | d£  dr/  df 


and  rewrite  the  stiffness  and  mass  matrices  as 


1 1 I 

[k]n  = J f f [B]  1 ID)  [B|  | I | djdrjdf  (3-14) 

-1  1 -1 


and 


1 1 1 

[m]  n = / f j [N]  1 Ip | [N 1 I J I df  dt/df  (3-15) 

1 1 -1 


for  Equation  3 1 3 to  be  satisfied  for  arbitrary  virtual  nodal  displacements  the  expression  in 
parentheses  must  vanish.  Invoking  this  condition  yields  the  equation  of  motion  for  element 
n which  can  be  written  using  Equations  3-14  and  3-15  as 


(3-16) 


IH 


I 


I quations  of  Motion  for  Linear  Viscoelasticity 

It  instead  of  Equation  3-1 3 we  specify  a linear  viscoelastic  relationship  between  stress 
and  strain,  considerably  more  complications  will  result,  but  the  complcxites  can  be  sur 
mounted  by  approximate  means.  The  result  is  that  the  head  injury  model  code  can  approx- 
imate the  time-dependent  relationship  between  stress  and  strain  exhibited  by  biological 
materials  when  such  a mechanism  is  justified. 

Since  considerable  attention  has  been  devoted  to  the  viscoelastic  nature  of  biological 
materials,  it  would  be  expected  that  a head  injury  model  must  possess  the  ability  to  account 
tor  viscoelastic  materials.  What  follows  is  a detailed  description  of  how  viscoelastic  material 
characterization  is  implemented  into  the  equations  of  motion.  1'his  work  is  largely  due  to 
Taylor  *s  . 

Viscoelastic  materials  are  often  called  “memory”  materials;  that  is,  the  current  state  of 
stress  in  the  material  is  determined  not  only  by  the  current  deformation,  but  also  by  all  past 
deformation  states.  Moreover,  the  memory  exhibits  a fading  phenomenon  in  that  past 
deformation  states  influence  the  current  stress  state  to  a lesser  degree  than  do  more  recent 
deformation  states.  I hese  characteristics  are  reflected  in  the  constitutive  model  for  visco- 
elasticity by  the  Sticltjcs  integral;  i.c., 


t(t){  [D(t-r)]  — |e(r)}dr 


(3-17) 


where  ^o(t)}  is  the  stress  vector  at  current  time  t,  je(r)}  is  the  strain  vector  for  0 < r < t,  and 
|l)(t  r ) 1 is  the  constitutive  matrix  composed  of  two  independent  relaxation  functions  for 
isotropic  materials.  This  constitutive  model  may  now  be  combined  with  the  eight-node 
finite  element  by  substituting  into  Equation  3 17  the  corresponding  strain-displacement 
relations  given  in  Equation  3-8.  Thus  we  have 


| o( t) | = /*[D(t-r)l  [ 15]  — 
or 


dr 


(3-18) 


JX.,  i -r  l 

K I laylor 
1966,  p.  21. 


An  approximate  method  tor  thermovisco-elastic  stress  analysis."  Nuclear  I ngmeering  ami  Design,  vol.  4. 
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Proceeding  as  before,  *his  expression  tor  Ju(t)}  is  substituted  into  the  virtual  work  equation 
tor  an  element  n (Kquation  3 12);  with  the  virtual  displacements  factored  out,  the  following 
is  obtained 


\gain  invoking  the  arbitrariness  of  the  virtual  displacements  in  the  above  equation  we  can 
write  the  equation  of  motion  for  a viscoelastic  element  as 


•11 

where  [k(t-r)|n  - f f f [13]  1 [D(t-r)]  [B]  |j|  d£  dr?  df 
1 -1  1 


fyi 


(3-20) 


and  [m] 11  is  as  defined  in  Equation  3-15.  Hence,  to  solve  the  general  linear  viscoelastic 
problem  (liquation  3-20)  utilizing  the  finite  element  method  requires  the  solution  of  simul- 
taneous linear  integral  equations.  Without  further  restriction  the  numerical  solution  of  in- 
tegral equations  of  this  form  requires  extensive  effort.  1 he  following  section  shows  a 
method  whereby  et  fort  can  be  greatly  reduced  without  undue  restriction  on  the  class  of 
problems  that  can  be  solved. 

Solution  of  Integral  Equations 

( onsidcr  a typical  integral  equation  governing  the  motion  of  a single-dcgree-of-frccdom 
s\  stem  with  viscoelastic  behavior 


t+  t 

mii  + / 

-oo 


(')u 


E(t  + t-T)  — 
<)t 


dr 


f(t+  t) 


(3-21) 
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It  the  body  is  in  equilibrium  until  time  zero,  Equation  3-21  becomes 

t+At  r)  u , 

mu  + E(t+  t)  u(o)  + f E(t+At-r)  — dr  = f(t+At)  (3-2*.) 

o ,)r 

It  the  relaxation  modulus  tor  the  material  is  characterized  by  a Prony  series,  a simplification 
to  the  solution  by  step  forward  time  integration  will  be  shown  to  result.  I herefore,  write 


E(t)  = E,.  + 


,1V  -t/X: 

“ EqC  * 

i=l 


(3-23) 


where  nv  is  the  number  of  terms  in  the  series. 

Substituting  liquation  3-23  into  l.quation  3 22  we  can  write 


nhi  + M0u(t+  t)  + 


-(t+At)/X- 

-'■e 


+ 


-(t+At)/X: 

K«e 


t+  ’ t 


t/X:  d u 
e 1 — dr 
'dr 


f(t+  t) 
(3-24) 


I his  equation  is  simplified  by  defining  the  last  term  on  the  left-hand  side  as  1;  i.e., 


I,(t+  t) 


-(t+  ' t)/X:  H 1 t/X 
**  ‘ ' - 


1 — dr 


( 3 25) 


and  can  be  expanded  for  subsequent  computational  advantage  as 


lj(t+  t)  e 


- t/X 


1 lj(t)  + Eje 


-(t+','t)/X: t+  1 t/X-  r)u 


< )t 


dr 


( 3 25a) 


l.quation  3-25a  is  a recursive  equation  for  step-forward  integration  of  Equation  3 24.  pro- 
vided the  integral  on  the  right-hand  side  ol  Equation  3 25a  can  be  solved  numerically,  lo 
do  this  it  is  assumed  that  the  displacement  u varies  linearly  over  each  time  increment.  In 
this  way 


9u  u(t+At)  - u(t)  u(t> 

f)T  At  At 

is  treated  as  a constant  over  the  interval  t<r<t+  t. 
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I he  fomulation  given  above  for  a singlc-degrcc-of-frccdom  system  holds  also  for  a sys- 
tem of  equations  for  any  viscoelastic  finite  element.  That  is 


22 


where  history  integrals  are  now  defined  by 


■(  1 j(t  + t)} 


-'•t/X: 


lj<t) 


+ [kj(  t)J 


' u(t)  ’ 
AV(t) 
i ■ ■w(t)> 


(3  2K) 


and  [k()|  and  [k,|  arc  viscoelastic  element  stiffness  matrices.  They  are  defined  in  the  same 
manner  as  element  stiffness  matrices  for  elastic  materials  (see  Equations  3-14  and  3 20). 


Normal  and  Reduced  Integration  of  l.lement  Stiffness  Matrix 


for  some  time  it  has  been  known  that  in  the  application  of  finite  element  methods  for 
nearly  incompressible  (i>  .=  0.5)  or  incompressible  materials  (u  = 0.5),  numerical  problems 
will  be  encountered  with  the  usual  displacement  theory.  To  discuss  how  such  materials  arc- 
dealt  with  in  the  ordinary  displacement  formulation  of  finite  element  analysis  it  is  neces- 
sary to  discuss  the  method  of  evaluating  the  integrals  defining  the  element  stiffness  matrices 
shown  m Kquations  3-14  or  3-20.  I he  procedure  is  equally  applicable  to  elastic  or  visco- 
elastic materials.  In  actual  practice  the  integration  over  the  element  volume  is  carried  out 
numerically  using  the  method  of  Gaussian  quadrature.  I he  numerical  form  of  liquation 
3-14  lor  an  element  n becomes 


nq  nq  nq 

[kln  = v v 2 Hj  [B]  1 11)1  [ B 1 IJI  (3-29) 

j=l  j=l  j=l 

In  this  triple  summation,  terms  on  the  right-hand  side  are  evaluated  at  specific  points 
within  the  element  known  as  Gauss  points.  The  II,  are  weighting  factors  and  their  values 
are  shown  in  I able  3-1 . Ordinarily  for  eight-node,  isoparametric  elements,  such  as  those 
employed  tor  the  skull  bone  in  this  investigation,  two-term  (nq  2)  quadrature  is  employed. 
This  implies  eight  Gaussian  points.  The  eight  Gauss  points  include  two  each  on  the  plus  and 
minus  sides  of  the  three  local  coordinate  axes  (see  figure  3 1 ).  However,  evaluating  the  ele 
ment  stiffness  matrix  for  nearly  incompressible  materials  is  another  matter. 

fried  '''  has  provided  mathematical  insight  into  why  the  normal  two  term  quadrature 
would  be  incorrect  for  nearly  incompressible  materials,  lie  shows  that  a reduced  integra- 
tion must  be  used  for  those  terms  involving  the  l.ame  constant  in  liquation  3-29  With  the 
eight-node  isoparametric  element,  the  correct  reduced  integration  is  one-point  Gaussian 


I I riol  "l;mnc  element  analysis  of  incompressible  material  by  residual  energy  balancing,"  International  Journal  of 
Solids  Structures,  vol  10,  1974,  pp  993  1002. 
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I able  3-1.  Abscissae  and  Weight  Coefficients  of  the 
Caussian  (Quadrature  f ormula 


1 

/ 

1 


I 

/ 

1 


1 

/ f(e,r?,f)de  dr?  df 

1 


nq  nq  nq 

“ “ “ ^ *iii,(alj  a2j  a3j^ 

J=1  J=1  J=1 


Number 

Abscissae, 

Weight  Coefficients, 

of  Terms, 

* a 

II 

nq 

1 

0 

2.0000  0000 

2 

0.57735  02691  89626 

1 .0000  0000 

3 

0.77459  66692  41482 

0.555  555  555 

0.00000  00000  00000 

0.888  888  888 

quadrature  at  the  element  centroid.  Naylor40  has  also  demonstrated  the  necessity  of  re- 
duced integration  in  several  finite  element  analyses  of  nearly  incompressible  materials.  This 
is  the  basic  technique  employed  for  brain  material  by  Shugar  and  Katona41  in  the  HIM 
code.  Though  it  can  be  employed  either  way,  the  HIM  code  currently  specifies  reduced  in- 
tegration tor  terms  involving  both  X and  (.  in  liquation  3-29.  No  instability  due  to  energy 
less  shear  modes  has  been  observed  when  using  reduced  integration  involving  all  terms  and 
when  the  elements  are  enclosed  or  confined. 

Solution  ot  the  liquations  of  Motion 

In  this  section  an  account  is  given  of  the  method  by  which  the  HIM  code  integrates  the 
dynamic  equations  of  motion  for  the  skull/brain  system  when  either  direct  or  indirect  im- 
pact loadings  are  specified.  I he  method  is  termed  a step-by-step  or  direct  integration  tech 
nique  and  is  suited  to  large  systems  of  equations  such  as  those  which  would  emanate  from  a 
three-dimensional  discretization  of  the  head. 


I)  I Naylor  Stresses  in  nearly  incompressible  materials  by  finite  elements  with  application  to  the  calculation  ot  ex- 
cess pore  pressures.'  intt  mational  Journal  lot  Numerical  Methods  in  I ngincering.  vol.  H.  1974.  pp  443-460. 

I \ Shugar  and  M (.  Katona.  "Development  ot  finite  clement  head  injury  model,  " Journal  of  I ngincering  Mechanics 
Division,  AS<  I vol  101,1  M3.  Jun  1975,  pp  223 -2  39 
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In  tlu-  realm  of  dynamic  finite  element  computer  codes,  much  has  been  written  concern- 
ing the  most  advantageous  method  of  numerically  integrating  the  equations  of  motion.  I he 
objective  in  this  investigation  was  to  canvass  the  available  integration  techniques  and  evalu- 
ate them  in  the  context  of  criteria  developed  for  large  three-dimensional  finite  element 
models  An  early  selection  of  an  operator  was  desirable  so  that  it  could  be  employed 
throughout  the  successive  model  development  and  its  behavior  under  various  conditions 
could  be  scrutinized  for  better  understanding  of  its  strengths  and  deficiencies. 

I he  criteria  for  selecting  a solution  scheme  are  much  the  same  for  all  finite  element 
codes4’  the  nature  of  the  problem  to  be  solved,  particularly  the  size  of  the  problem  and 
the  load-time  characteristics. 

I he  first  step  is  to  choose  between  a modal  superposition  approach  or  a direct  integra- 
tion technique  I he  former  approach  is  not  feasible  tor  large  problems  unless  only  a few  of 
the  lower  modes  are  desired,  which  implies  that  the  loading  is  sluggish  and  only  excites  the 
lower  modes  * ' 44  In  the  case  of  the  head  injury  model  the  problem  is  very  large.  Within 
the  realm  of  direct  impact  loading,  short  load  durations  arc  the  rule  rather  than  the  excep- 
tion and  one  can  reasonably  expect  to  encounter  the  wide  spectrum  of  frequencies.  Thus, 
the  loading  can  excite  substantially  a number  of  the  higher  modes,  f urthermore,  the  mode 
superposition  approach  presumes  linearity.  Therefore,  because  the  system  is  large,  because 
the  loads  will  possess  short  durations,  and  because  an  extension  to  nonlinear  behavior  is  en 
visioned,  the  modal  superposition  approach  was  discarded  in  favor  of  a direct  integration 
technique  I lus  decision  also  related  to  the  selection  of  a finite  element  computer  code  lor 
modeling  the  head  \ modal  analysis  capability  was  not  considered  relevant  in  the  code 
selection  process 

I he  choice  is  then  reduced  to  that  of  either  an  implicit  or  an  explicit  integration  tech- 
nique I he  explicit  method  has  the  advantage  of  being  extremely  fast  computationally  for 
each  time  step  because  no  triangularization  is  involved.  However,  the  time  step  size  is  re- 
stricted b\  the  ( ourant  stability  criterion.  I hat  is,  the  time  step  must  be  less  than  the 
smallest  sonic  travel  time  across  any  element.  Unless  some  means  is  provided  within  the 
code  to  control  the  time  step  size  during  execution,  the  user  runs  the  risk  of  an  unstable 
solution,  for  codes  dealing  with  large  three-dimensional  discretizations  intended  for  re- 
peated use  in  parametric  studies,  the  possibility  of  instability  is  to  be  avoided.  I T en  w ith 


4 ’ 

" l<  S Dunham,  K I Nickell.  .iml  D < Stickler  ‘Integration  operators  lor  transient  structural  response Comp  tiers 
anti  Structures,  vol  2,  no.  1/2,  Ich  1*272  pp  1 16 

* S<  W ('.lough  and  J Pen/ien.  Dynamics  ot  structures.  New  York,  McCiraw  Hill  Hook  Company,  1975  p 271 
44 

l<  I)  < nok  ( oncepts  an. I applications  ot  finite  element  analysis  New  York.  John  Wiley  anil  Sons,  Inc  , 1974.  pp 
252-253. 
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permissible  time  steps  an  exorbitant  number  would  be  required  to  span  typical  load  dura- 
tions These  considerations,  which  arc  primarily  practical  in  nature,  reduce  the  selection  to 
one  ot  the  implicit  integration  operator  techniques.  Some  of  these  techniques  arc  uncondi- 
tionally stable  anil  pros  ide  control  over  computer  costs  in  allowing  a choice  of  time  step 
si/e  independent  ol  the  Courant  stability  criterion. 

Among  the  implicit  schemes,  the  Newmark  (7  1/2.  1/4)  operator45,  the  Wilson 

averaging  operator46,  and  the  Houbolt  operator47  are  more  prevalent.  I he  accuracy  of 
these  methods  was  studied  to  decide  which  would  be  more  suitable  for  the  head  injury 
model 

Inherent  in  all  implicit  schemes  is  the  tendency  to  artificially  attenuate  the  response  and 
tn  artificial!)  elongate  the  period  and  shift  the  phase  of  the  response  during  integration  of 
the  equations  of  motion. 

Response  attenuation  for  each  operator  was  investigated  with  a simple  1-degree-ol  - 
treedom  linear  oscillator.  This  procedure  was  first  used  to  study  stability  of  finite  differ- 
ence solutions4'4  and  has  been  used  contemporaneously  by  others49  in  evaluating  the  rela 
tive  merits  of  integration  operators.  Only  a heuristic  account  of  the  evaluation  procedure  is 
presented  in  this  paper.  More  complete  information  is  contained  in  the  references  cited. 

I he  free  response  of  a single-degrec-of  freedom  linear  oscillator  is  governed  by  the  dif- 
ferential equation  of  motion 


••  2 

u + co  u 


t) 


(3-30) 


where  u is  the  displacement  and  us  is  the  natural  circular  frequency.  The  procedure  begins 
b>  substituting  the  integration  operator  expressions  presented  in  fable  3-2  into  1 quation 
3 3o  .11  time  tn+  | In  this  wa\  a numerical  form  of  the  equation  (also  shown  in  I able  3-2) 
tor  each  operator  is  derived,  and  it  w ill  be  necessary  to  find  the  expression  for  a typical 


Nathan  VI  Neumark.  \ method  of  computation  for  structural  dynamics,"  Journal  of  the  I ngineermg  Mechanics 
Division.  AS<  I . I ,MJ.  Jul  195‘> 

niversity  of  < aliforrna  SI  s\1  Report  I A computer  program  for  the  dynamic  stress  analysis  of  underground 
structures.  Iiy  I I Wilson.  Berkeley,  < alit.Jan  |V6K. 

|ohn  c llouholt  \ recurrence  matrix  solution  for  the  dynamic  response  of  elastic  aircraft, " Journal  of  Aero  Science, 
vol  I 7.  1 mso.  p 540 


I*  I)  I as  and  K I)  Kichtmeyer  Surve  ' of  the  stability  of  linear  finite  difference  equations,”  Communications  on 
I’vm  and  Applied  Mathematics,  vol  9.  no  2,  Ma\  1956 
49 

Robert  I Nickcll  On  the  stability  of  'pproximation  operators  m problems  of  structural  dynamics,"  International 
Journal  ot  Solids  and  Structures,  vol.  7.  1 ‘>7  I . pp  301  3 I'/ 
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I able  3-2.  Integration  Operator  Description 
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response  un+  j for  each.  After  the  particular  operator  under  evaluation  has  been  introduced, 
the  resulting  numerical  equation  of  motion  is  then  east  into  a difference  equation  of  motion 
by  using  difference  approximations  to  the  derivatives.  The  oscillating  component  of  the 
solution  (in  polar  form)  of  this  difference  equation  is, 


un+  1 


,n  + 1 


A 1 c 


,iu(n+  1 ) 


A 2e'ia<n+ 1 ) 


(331) 


I he  terms  p .uni  u are  the  polar  modulus  and  polar  phase,  respectively,  and  the  constants 
A I anil  \i  are  determined  from  initial  conditions. 

f rom  Kpiation  3 31  it  is  seen  that  the  response  will  grow  without  limit  and  therefore 
he  unstable  it  /;  > 1 However,  it  p < 1 the  response  will  be  attenuated  continuously.  I.ach 
ot  the  operators  will  Yield  a different  function  lor  the  modulus  p in  terms  of  a parameter 
0 which  in  turn  is  the  product  of  frequency  oj  and  time  step  si/e  t Although  attenuation 
occurs  for  p 1 , this  condition  is  preferable  to  instability  and  is  the  requirement  for  a 
stable  integration  operator.  1 he  results  of  this  investigation  are  presented  in  f igure  3-2, 
where  the  moduli  p are  plotted  as  a function  of  <;>.  It  is  seen  that  all  three  operators  are  un- 
conditionally stable;  i.e.,  p I for  any  time  step  size  t.  further,  it  is  seen  that  the  New 
mark  operator  possesses  no  attenuation  and,  on  that  basis,  is  the  choice  among  the  three 
operators  investigated,  t his  conclusion  is  corroborated  in  Reference  42  which  states  further 
that  the  New  mark  operator  is  also  a good  choice  for  nonlinear  applications. 

Other  considerations  could  be  important  criteria,  depending  on  the  investigator’s  objec- 
tives. l or  example,  computer  time  (number  of  numerical  operations)  is  important  although 
it  is  believed  that  the  three  operators  studied  are  economically  equal,  f urthermore,  an 
evaluation  based  on  a single  dcgrec-of-frecdom  system  may  be  inconclusive  inasmuch  as 
subsequent  application  of  the  operator  is  with  highly  complex  systems,  fo  the  writer’s 
knowledge  there  are  no  reliable  evaluation  techniques  which  account  tor  system  complex- 
ity. I he  New  mark  operator  (d  1/4,7  1/2)  will  evidence  numerical  damping  in  systems 

more  complex  than  one-degrec-of ’-freedom  (i.e.,  general  systems). 

Instead  of  the  entire  complex  system  of  the  skull,  the  mathematical  focus  here  will  be 
on  the  integration  of  the  equations  associated  with  two  typical  arbitrary  mass  points  (or 
nodes)  within  the  global  system,  as  shown  in  f igure  3-3.  from  this,  the  salient  features  ot 
the  solution  process  can  be  illustrated  without  direct  reference  to  the  remainder  of  the 
system. 

DIKf  C I IMPAC  I I he  governing  equations  of  motion  for  arbitrary  mass  points  i and 
i+1  at  time  step  Sr  1 when  subjected  to  specified  direct  impact  forces  P(r)  are 
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I' igurc  3 3 C.lnlial  Displacement  of  Skull  liram  Model  Relative  to  Inertial  Reference  Sv stem 


I he  mass  terms  nij  and  mj+]  represent  lumped  masses  at  the  two  nodes  in  a diagonal  mass 
matrix  1 he  global  stiffness  terifis  j are  obtained  in  the  direct  stiffness  method14  by  an 
appropriate  summation  of  element  stiffness  terms  from  Equation  3-29.  That  is,  j 
-kj|.  1 he  displacement  vector  and  acceleration  vector  composed  of  Uj  and  LJ ■ terms,  re- 
spectively, are  unknown  at  time  step  S+l.  By  definition,  they  are  referred  to  the  inertial 
reterenee  system.  By  choosing  sufficiently  small  time  increment  t,  it  is  possible  to  predict 
the  acceleration  turns  at  step  S+l  in  terms  of  known  quantities  from  previous  time  steps 
and  the  unknown  L- 1 1 1 he  particular  prediction  algorithm  or  integration  operator  em- 

ployed, as  previously  described,  is  the  Newmark  d method. 


|:.(S+1) 

i 


ir.(S+  1 ) 


|V(S) 


At 


uVsV'tI 2 


) 


(3-33) 


Substituting  Equation  3-33  into  Equation  3-32  and  combining  the  unknown  Uj  terms  on  the 
left  hand  side  anil  taking  the  known  terms  to  the  right-hand  side,  the  following  equation 
results. 
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I he  modified  force  vector  on  the  right-hand  side  has  been  updated  with  the  addition  of  in- 
ertia forces.  I he  modified  stiffness  matrix  on  the  left-hand  side  never  requires  updating  for 

linear  analyses  with  constant  time  increments.  It,  therefore,  needs  to  be  inverted  only  once 
(at  the  beginning,  prior  to  stepping  through  time).  The  left-  and  right-hand  sides  of  Equa- 
tion 3-34  are  then  premultiplied  by  the  inverse  to  obtain  the  solution  for  unknown  dis- 
placements at  step  S+l . I lowever,  inverting  and  multiplication  by  the  inverse  as  described 
here  is  only  symbolic  In  actual  practice,  the  matrices  arc  too  large  for  inversion  to  be 
economical,  and  instead  the  procedure  of  Gaussian  elimination  is  employed  to  solve  Equa- 
tion 3 34  for  the  unknown  displacements. 
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INIJIRI  (.  1 I Ml'  At  I lo  simulate  indirect  impact  loads  the  displacement-time  history 
resulting  trom  impact  is  specified  instead  of  force  histories  I or  example,  it  the  forces  in 
I igurc  3-3  are  assumed  zero  and  the  displacement  Uj+  ,(t>  of  the  i+  1 degree  of  freedom  is 
assumed  known.  I quation  3 32  mav  he  written  as. 
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Using  the  Newmark  operator  by  substituting  liquation  3 3 3 into  liquation  3-35  we  obtain, 
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Isolating  the  ith  equation  ot  I quations  3-36  anil  rewriting  it  by  placing  all  the  known  terms 
on  the  right-hand  side,  we  obtain 
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Since  l j+j  is  known,  the  equation  associated  with  it  (the  second  ot  Equations  3-36)  can 
be  eliminated  from  the  system  and  Equation  3-36  may  be  rewritten  as 
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Again,  in  practice,  liquation  3 38  is  solved  for  the  unknown  Uj  by  Gaussian  elimination 
rather  than  by  actual  inversion  of  the  square  matrix  of  the  left-hand  side. 

When  p is  assigned  a value  of  1/6  in  Equation  3-38,  the  integration  of  the  equations  of 
motion  will  be  only  conditionally  stable,  and  the  linear  acceleration  method  results.  A! 
though  the  parameter  p is  optional  in  the  HIM  code,  a value  of  1/4  has  been  exclusively  used 
in  the  development  of  the  head  injury  model.  In  this  case,  the  method  is  unconditionally 
stable  and  is  termed  the  average  or  constant  acceleration  method.  In  general  neither  method 
is  considered  more  accurate  than  the  other.  With  the  use  of  p 1/4  the  analvst  does  not  run 
the  risk  of  having  the  solution  grow  without  bound  for  an  arbitrary  t 

The  direct  integration  technique  also  provides  flexibility  in  handling  a variety  of  arbi- 
trary impact  loads  whose  frequency  content  may  be  quite  varied  and  extensive  It  is  this 
frequency  content  that  determines  which  skull/brain  modes  need  to  be  integrated  accurately 
and,  hence,  what  si/e  of  time  step  is  required. 
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At  the  beginning  of  the  head  injury  model  development  program,  no  doeuniented  evi 
ilcncc  existed  to  indicate  that  the  finite  element  method  was  a suitable  analysis  tool  only 
the  suspicion  that  it  might  contribute  to  the  study  of'  head  injury  due  to  its  innate  ability  to 
structurally  model  solid  forms  with  complicated  shapes  and  constitutive  properties. 

Manx  uncertainties  about  its  applicability  to  head  injury  modeling  existed,  not  the 
least  of  which  was  due  to  the  known  semisolid  or  fluidlike  makeup  of  brain  matter. 

Consistent  with  a cautious  philosophy,  a series  of  one-  and  two-dimensional  models  were 
constructed  and  numerically  studied  prior  to  attempting  a three-dimensional  model.  These 
studies  are  reported  in  this  section,  however  some  readers  may  wish  to  skip  directly 
to  the  three-dimensional  studies  presented  in  section  5. 


One  Dimensional  Models 

Since  it  was  originally  envisioned  that  an  assemblage  of  eight-node  brick  finite  elements 
would  constitute  the  final  three-dimensional  model  configuration,  the  eight-node  brick  ele- 
ment was  chosen  to  comprise  the  basic  one-dimensional  model  used  in  this  study. 

With  the  proper  choice  of  bulk  and  shear  moduli,  a Poisson  ratio  of  /.ere;  can  be  obtained 
and,  in  effect,  produces  a one-dimensional  model  from  a prismatic  stack  of  three- 
dimensional  elements  as  shown  in  f igure  4-1  Alternatively,  boundary  conditions  can  be 
employed  to  establish  the  same  effect  and  retain  the  freedom  to  specify  arbitrary  bulk  and 
shear  moduli  values.  I he  first  of  these  methods  was  employed  in  this  study  to  facilitate  an 
evaluation  of  the  eight-node  brick  by  comparison  with  classical  one-dimensional  wave  propa- 
gation theory.  I hc  second  method  was  used  in  the  latter  part  of  the  study  where  the  use  of 
realistic  skull/brain  material  properties  was  necessary.  Specifically,  boundary  conditions 
were  assigned  to  restrain,  completely,  lateral  deformations  (transverse  to  direction  of  load- 
ing) in  the  one-dimensional  stack  of  elements.  This  was  necessary  to  prevent  physical  in- 
stability induced  by  the  zero  value  assigned  to  the  shear  modulus  for  brain  material.  I he 
stack  of  three-dimensional  eight-node  brick  elements  was  subjected  to  short  stress  pulses 
typically  found  in  stress  wave  problems.  It  is  important  to  point  out  that  pulse  durations 
typical  of  head  injury  are  very  often  larger,  producing  wavelengths  much  in  excess  of  cranial 
dimensions.  As  a result,  the  head  injury  problem  belongs,  most  often,  in  the  vibration  do- 
main, and  the  code’s  ability  to  propagate  stress  waves  accurately  is  not  of  primary  signifi- 
cance. Nevertheless  many  engineers  heretofore  have  investigated  the  dynamics  of  head 
injury  from  the  stress  wave  standpoint;  their  main  concern  being  the  potentially  injurious 
effect  of  discrete  tensile  waves  in  brain  matter  if  and  when  they  occur  (again,  they  owe  their 
existence  to  i -datively  small  pulse  durations),  because  of  the  influencing  nature  of  these 
past  investigations  and  because  it  is  conceivable  that  small  pulse  durations  may  be  found  to 
exist  among  causative  head  injury  loads,  it  is  desirable  to  know  how  well  the  intended  for- 
mulation performs  in  such  instances. 

An  individualized  summary  of  each  model  used  in  the  study  follows.  The  information 
includes  geometry,  material  properties,  input  loads,  acoustic  velocity  and  critical  time  step 
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values  as  well  as  the  computed  stress  wave  results.  It  is  believed  that  a one-dimensional 
model  does  not  adequately  represent  the  cranium;  thus,  no  conclusions  contained  in  this 
section  are  intended  to  be  applicable  to  head  injury. 

Fguilibrium  Model.  I he  primary  objective  of  the  equilibrium  model  was  to  verity  that 
the  I I AT  code  and  particularly  the  eight-node  brick  were  operating  correctly.  An  cquili 
brium  check  in  finite  clement  work  is  usually  the  means  used  to  make  such  a verification.  I! 
a 1 0 psi  step  pulse  load  is  input  as  shown  in  Figure  4-2,  the  computed  element  normal  stress 
<;vv  should  take  on  the  value  of  1 .0  psi.  Figure  4-3  shows  this  stress  plotted  against  time 
lor  three  different  time  steps  and  for  each  of  the  four  elements  of  the  equilibrium  model. 
These  results,  especially  those  for  element  number  1,  demonstrate  that  the  code  is  indeed 
operating  correctly.  Oscillation  of  the  computed  stress  about  the  theoretical  value  is  to  be 
expected  with  any  numerical  integration  scheme  whether  it  is  implicit  or  explicit. 

Because  the  stress  wave  is  reflected  (with  the  same  sign)  from  the  fixed  end  face  of  the 
model,  stresses  are  correctly  showing  a tendency  in  elements  2 through  4 to  increase  to  a 
value  of  -2.0  psi. 

More  importantly,  the  effect  of  time  step  size  on  accuracy  is  illustrated.  ( liven  sufti 
cient  time  ( 1 second)  after  the  initiation  of  loading,  data  for  the  larger  time  step  t 0.8 
second  (approximately  twice  the  critical  value)  appears  as  accurate  as  data  for  smaller  time- 
steps.  In  fact,  the  greater  numerical  damping  associated  with  this  larger  time  step  actually 
appears  beneficial  in  damping  the  oscillation. 

Reflection/ Refraction  Model.  Because  it  is  of  interest  to  know  how  well  the  11  AP  code- 
handles  layered  media,  a simple  model  composed  of  two  materials  was  numerically  tested 
prior  to  using  more  realistic  values  for  material  and  material  layer  thickness.  The  model  is 
shown  in  Figure  4-4. 

According  to  one-dimensional  stress  wave  theory,  a wave  incident  upon  a material  inter- 
face where  the  characteristic  impedance  changes  abruptly  will  be  both  reflected  and  re- 
tracted. This  theoretical  behavior  applied  to  the  present  model  was  computed  and  is  show  n 
as  the  dashed  line  response  in  Figure  4-5. 

It  can  be  observed  that  the  computed  solution  is  responsive  to  the  reflection  and  refrac- 
tion of  waves.  In  Material  1 , the  average  of  the  three  greatest  computed  values  for  times  of 
2 4,  3.2  and  4.0  seconds  appears  to  oscillate  about  -1.3  psi,  which  is  close  to  the  magnitude 
of  the  returning  (reflected)  wave.  Since  Material  2 involves  a steep  wave  front,  agreement  is 
not  good,  but  still  the  average  of  the  computed  refracted  wave  increases  with  time  toward 
the  theoretical  value. 

Numerical  Damping  Model.  This  model  was  intended  for  investigating  the  numerical 
damping  which  will  occur  with  the  Ncwmark  jj  method  when  it  is  employed  in  systems  more 
complex  than  a sing!e-dcgree-of-frccdom  system.  The  model  is  shown  in  F igure  4-6. 

I he  results  of  this  model  clearly  show  the  numerical  damping  effect  of  large  time  steps 
on  the  stress  wave  solution.  Figures  4-7,  4-8,  and  4-9  each  illustrate  the  stress  wave  solution 
at  various  times  as  computed  by  the  Fl-'.AP  code  for  three  different  time  step  sizes ; they  are 
time  steps  equal  to  the  critical  value,  twice  the  critical  value,  and  four  times  the  critical 
value. 
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l-igure  4-4.  Kefleet ion/ Kefraetion  Model  Description 


\S  .1  result  ot  this  stuilv,  .in  iile .1  can  he  gained  as  to  the  relame  merits  ot  the  cost  time 
step  si/e  tradeotl  ( .cncrally  . the  smaller  the  time  step,  the  more  accurate  and  costly  the 
,mah  sis  u til  he  I Ins  is  characteristic  ot  am  implicit  scheme,  but  more  importantly . the 
tradeotl  betw ci  n cost  and  .u  curai \ can  he  controlled  On  the  other  hand,  an  explicit  inte- 
gration si  heme  would  he  unstable  lor  am  time  step  si/e  greater  t han  the  critical  \ aluc  and 
thus  no  such  control  would  exist 

l hese  results  suggest  that  a time  step  si/e  greater  than  critical  hut  less  than  tw  ice  critical 
max  Ik  accept  able  in  am  subsequent  analyses  invoking  discrete  waves  ( main!)  it  appears 
that  too  much  damping  is  present  with  a time  step  size  equal  to  lour  t lines  the  critical  value 

I wo  additional  observations  can  he  made  I irst,  numerical  damping  retards  the  vv ave 
propagation  speed,  i his  can  he  seen,  tor  example,  by  comparing  the  magnitude  ot  the  lead 
ing  edge  ot  the  wave  at  the  3 inch  station  along  the  column  tor  the  three  dittcrcnt  time 
steps  I his  value  is  highest  in  I igurc  4-4  and  lowest  in  I igurc  4-6  Secondly,  the  tensili 
stresses  in  the  trailing  edge  ot  the  waves  are  due  to  numeric.il  noise  and  become  more  pm 
nounced  with  larger  time  step  sizes 

I av  en  d Media  Mod<  I \ triangular  pulse  identical  to  the  one  used  in  the  prev ious 
model  was  appiieJTo  a lav  end  media  model  as  shown  in  I igurc  4 1<>  I he  only  ditlerence 
between  the  two  models  is  t hi  layered  bone  caps  at  both  ends  ol  the  stack  ol  brain  elements 
that  are  added  to  the  present  model  I herelore,  a direct  comparison  ol  I igurc  4-7  and  the 
results  for  this  model  ( I igurc  4 11)  illustrates  the  influence  ot  the  hone  caps  I hi  hit  i I m 
both  computer  runs  were  bused  on  the  same  time  step  si/e  with  relatively  little  ch.mgi  to 
brain  element  discretization. 

Overall  the  elastic  hone  caps  appear  to  reduce  slightly  the  subsequent  peak  stresses  in 
the  brain  matter  and  subject  the  brain  elements  to  an  earlier  pulse  Otherwise,  the  el  led 
nl  the  elastic  bone  characterization  on  transmitted  stress  waves  is  very  minimal  in  this  study 
I lie  slight  stress  reduction  could  be  attributed  to  numerical  damping  that  vv  as  not  present  in 
I igurc  4-7  but  may  exist  in  the  present  model  bei.iuse  the  time  step  is  much  larger  <bv 
about  7 times)  than  a computed  critical  value  lor  the  bone  elements 

Summarv  ol  t )iu  I )i im  i isi< . . . , ! Model  Study  Kcsnlts  familiarity  with  the  c h .tract eristics 
ol  input  data  format  and,  more  importantly . vv  it h computed  output  data  using  the  I I \P 
code  vv  as  gamed.  I he  st ml v was  dccidedlv  oriented  to  the  vv  av  e me.  names  realm  and  the 
implicit  technique  ol  integration  evaluated  within  this  cm  ironment  In  am  future  analy  ses 
with  the  II  \ I’  code  that  in. iv  require  small  time  steps,  it  is  expci ted  t hat  accurate  results 
using  near  critical  time  steps  can  be  achieved  \lso,  in  analyses  involving  problems  w here 
wavelengths  are  on  the  order  of  i ramal  dimensions  and  beyond . t is  expei ted  that  accurate 
results  can  again  tie  achieved  with  a corresponding  savings  in  cost  through  the  use  ot  turn 
steps  substantial  greater  than  critical  values  < ontrol  over  the  decision  between  cost  and 
accuracy  is  provided. 

Specific  conclusions  resulting  from  this  study  follow 

I I quihbrium  checks  on  a nominal  one  dimensional  homogeneous  model  with  step 
loading  showed  that  the  I I \ I * code  is  capable  ol  matching  trends  of  classical  theory 
solutions 
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2 Reflcction/rcfraction  tests  on  a nominal  one-dimensional  dissimilar  material  model 
with  step  loading  showed  that,  on  the  average,  the  I I.AP  code  captured  well  the  reflected 
wave  magnitude.  Results  for  the  refracted  wave  were  less  accurate  due  to  the  severity  of  the 
wavefront  but  the  tendency  to  capture  it  was  evident 

3 Investigation  of  the  effect  of  time  step  si/e  on  a one-dimensional  brain  model  with  a 
triangular  pulse  loading  of  short  duration  showed  that  critical  time  step  sizes  yield  results  in 
excellent  agreement  with  theoretical  prediction  However,  as  the  prescribed  time  step  si/c 
increases,  significant  numerical  damping  increases  and  adversely  affects  the  solution.  \Iso, 
the  use  of  larger  time  steps  resulted  in  a retardation  of  w ave  propagation  speed  and  an  in 
crease  in  noise  in  the  trailing  portion  of  the  waves 

4 I'he  addition  of  elastic  layered  bone  caps  to  the  one-dimensional  brain  model  had 
negligible  eff  ect  on  the  stress  wave  response  in  the  brain 

\xixs  mmctric  Model 

I he  objective  of  this  model  was  to  determine  whether  or  not  the  finite  element  method 
was  able  to  function  as  a head  injure'  modeling  technique.  Manv  investigators  have  used, 
and  arc  still  using,  a.xisymmetric  analytical  models  for  head  injury  modeling,  and  it  was  felt 
that  if  the  finite  element  method  could  duplicate  some  results  previously  obtained  from  the 
models  it  too  could  be  used  as  a viable  analysis  tool  and  be,  in  principle,  easily'  extended  to 
more  realistic  head  injury  modeling,  liven  though  the  geometry  is  simplified  in  the  axisvm 
metric  form,  tin-  ability  of  the  finite  element  method  to  handle  a very  soft  material  en- 
capsulated by  a very  stiff  (albeit  not  rigid ) material  w as  of  primary  concern  and  required 
demonstration. 

A finite  difference  study  by  Merchant  and  (irispino'"  was  used  in  the  comparison  w ith 
the  finite  element  analysis.  I hesc  investigators,  in  turn,  found  good  agreement  between 
their  model  and  the  closed  form  work  of  I ngin1'  . 

I he  design  of  the  a.xisymmetric  finite  element  model  was  somewhat  affected  bv  con 
dieting  requirements  on  one  hand,  the  objective  was  to  compare  w ith  a f inite  difference 
solution;  on  the  other,  it  was  desired  to  demonstrate  the  model  as  a valid  head  injurs  study 
tool  in  its  own  right.  I he  finite  difference  study  contained  a two -layered  skull.  I he  present 
writer  believed  that  a three  layered  skidl  representing  the  table/d iploe" complex  was  manda- 
tors , based  upon  the  work  of  Hards  and  Marcal1*,  and  therefore  implemented  the  three- 
layer  simulation.  Despite  the  dif  ference,  the  models  are  beliesed  to  be  of  sufficient  similar- 
ity lor  comparative  purposes  since  the  diameters,  degree  of  discretization,  and  the  material 
properties  were  the  same  for  both  models. 


I mvtrsiiyol  W ishmgfon.  Department  of  Wch.innal  I nginctrm^  \ dynamic  analysis  I an  clastic  model  of  the  human 
he.nl.  I » y Howard  ( Merchant  and  Anthoin  I < rispino  Seattle.  Wash,  l*/  '.1 

* All  I I ngm  I he  axisymnictric  response  id  i fluid  filled  spherical  shell  to  a loc  al  radial  impulse  a model  for  head 
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I he  finite  dement  model  can  be  seen  in  f igure  4-12.  I lie  brain  mesh  was  obtained  by 
means  of  a I aplacian  generation  scheme  where  ail  unspecified  nodal  points  are  computed 
iteratively  until  they  are  uniformly  spaced  and  consistent  with  specified  node  points  tor  the 
skull  I his  technique  is  discussed  in  more  detail  in  section  5. 

Material  properties  are  summarized  in  fable  4-1 . Different  dimensional  systems  were 
employed,  but  all  v alues  are  equivalent  I he  tabulated  mechanical  and  physical  properties 
for  hard  bone  and  brain  are  traceable  to  ( ioldsmith5’  I lowever,  for  diploe*,  material  proper- 
ties were  computed  as  suggested  by  Melvin,  et  al.53.  Namely,  bulk  and  shear  muduli  were 
calculated  by  multiplying  their  respective  values  for  hard  bone  by  the  weight  density  ratio 
between  diploe’ anil  hard  bone.  In  this  case  the  ratio  was  one-tenth  and  the  resulting  bulk 
and  shear  moduli  for  diplo’c" are  1 333  x 10^  psi  and  0.80  x HP  psi , respectively  Note  that 
with  this  ratio  technique  the  elastic  sonic  velocity  is  the  same  for  both  hard  bone  and 
diploe’. 

Prior  to  conducting  dynamic  studies  with  this  model,  a static  axisymmetric  load  was  ap- 
plied radially  at  one  pole  of  the  sphere  while  restraining  the  opposite  pole  against  all  dis- 
placement. A computed  hydrostatic  pressure  state  was  expected  in  this  simulation  but  did 
not  occur  despite  the  mechanical  properties  for  the  internal  elements  being  specified  as 
those  for  a fluid  (see  I able  4-1 ).  What  was  computed  was  a highly  variable  and  random 
pressure  distribution  (predicted  magnitudes  varied  from  +6  to  -194  psi)  within  the  sphere. 

I his  apparent  anomaly  was  really  a misinterpretation  of  the  simulation  because  the  outer 
brain  elements  are  attached  everywhere  to  the  inner  skull  elements,  in  contrast  to  the  actual 
behavior  of  fluid  that  would  not  attach  itself  to  the  skull  elements:  Thus,  the  computation 
w as  correct  but  the  initial  interpretation  of  the  simulation  was  incorrect  since,  in  this  case, 
the  desire  was  to  create  a hydrostatic  pressure.  Means  were  sought  to  eliminate  the  inherent 
and  troublesome  displacement  continuity  between  the  skull  and  brain  finite  elements  The 
difficulty  was  that  bending  and  shear  modes  were  being  transferred  from  the  shell  elements 
into  the  interior  elements  because  of  nodal  point  compatibility  at  the  shell/interior  inter- 
face. It  was  believed  that  the  model  must  include  a mechanism  to  simulate  the  apparent 
discontinuity  between  the  skull  and  brain  that  arises  due  to  the  presence  of  the  subarach- 
noid space.  I here  fore,  the  model  must  either  relinquish  compatibility  between  the  nodes  of 
the  shell  and  the  nodes  of  the  core  or  prov  ide  a substitute  mechanism.  Detaching  or  co- 
alescing nodes  at  a common  interface  in  thejfinitc  element  method  can  become  difficult  and 
expensive  where  a large  number  of  nodes  art  concerned,  t herefore,  an  alternative  approach 
was  taken  which  resulted  in  the  development  of  a brain  element. 
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I able  4 1 . Mechanical  and  Physical  Properties  lor  Axisymmetric  Models 


Property 

Kngin  and 
( ioldsniith^ 

Crispino  and 
Merchant3 

<1.1 

a.  Bone  (inn 

er  and  outer  layer) 

P 

0.0772  lbm/in.3 
(2.14  g/ent®) 

2.14  g/cm  ® 

-4  lbf-s2 

2.0  x 10 

in  4 

k 

(1.333  x 106  psi) 

9.23  x H)4  bars 
(1.333  x 106  psi) 

1.3  33  x 106  psi 

(i 

(0.8  x 106  psi) 

5 .5  3 x 1 04  bars 
(0.80  x 106  psi) 

0.8  x 106  psi 

<1 

(109,500  in./s 

278.5  cm/ms 
( 109,500  in./s) 

109,500  in./s 

1 >* 

2 x 1 06  psi 

(2  x H)6  psi) 

2 x 1(/’  psi 

•v 

0.25 

(0.25) 

0.25 

b.  Brain 

P 

0.0362  lbm/in.® 
( 1 .0  g/cm®) 

1 .0  gm/em® 

4 lbt-s 
0.937  x 10 

in  4 

K 

(0.305  x 106  psi) 

1.96  x 104  bars 
(0.084  x 106  psi) 

0.305  x 106  psi 

(, 

(0) 

0 

0 

<d 

57,100  in./s 

1 44  cm/ms 
56,700  in./s 

57,100  in./s 

(0) 

(0) 

0 

7 

(0.5) 

(0.5) 

0.5 

’I  2 x 106  psi  is  lor  hard  bone;  for  diploe’,  I 2 x 10®  psi  (assumed). 
'Values  in  parentheses  are  converted  front  authors’  values. 


I 


I he  material  eharaeteristies  of  brain  matter  have  been  reported  and  are  generally  ae 
eepted  as  exhibiting  a nearly  incompressible,  constant  bulk  modulus,  and  a very  small  time- 
dependent  shear  modulus. 

In  Reference  54  the  use  of  a very  compressible,  “effective,”  bulk  modulus  for  brain 
material  is  discussed  In  many  eases  of  finite  element  modeling,  the  use  of  effective  material 
property  constants  is  justifiable  because  the  actual  values  are  either  unknown  or  because, 
through  their  specification,  the  solution  becomes  intractable.  If  brain  injure  is  to  be 
modeled,  the  correct  and  unique  value  relating  a state  of  strain  to  a state  of  stress  in  brain 
material  is  dcsireablc.  I his  relationship  can  be  provided  by  a high  (liquid)  bulk  modulus 
( * 305,000  psi).  An  “effective”  bulk  modulus  which  results  from  observed  distensibility  of 
the  entire  central  nervous  system,  or  even  the  cranial  compartment  itself  may  not  necessar- 
ily lie  required.  I lie  distensibility  mechanism  can  be  provided,  for  example,  by  simulating 
a compressible  subarachnoid  space. 

t he  characteristics  of  an  ideal  fluid  imply  the  presence  of  hydrostatic  stress  proportional 
to  volume  change  and  the  absence  ol  shear  stresses.  Front  a classical  mechanics  viewpoint, 
a fluid  can  be  characterized  very  simply  by  specifying  the  appropriate  bulk  modulus  and 
setting  the  shear  modulus  to  zero  for  those  regions  representing  the  fluid.  However,  when 
employing  the  finite  element  technique,  further  precautions  must  be  taken  to  avoid  erro- 
neous results.  Specifically,  the  standard  use  of  high  order  numerical  integration  techniques 
to  obtain  the  clement  stiffness  matrix  must  be  replaced  by  a reduced  numerical  integration 
(in  this  ease,  a first-order  integration)  capable  of  determining  the  volume  change  of  the  ele- 
ment independent  of  the  assumed  interpolation  functions  for  the  displacements. 

I he  reason  for  this  reduced  integration  is  to  insure  that  the  strain  energy  of  the  element 
responds  only  to  deformation  modes  that  are  associated  w ith  volume  change  rather  than 
deformation  modes  such  as  shear  and  bending.  Not  only  did  this  approach  work  well  but 
also  proved  a very  inexpensive  coding  modification.  From  this,  it  was  learned  how  a brain 
element  could  be  provided.  In  summary  the  brain  element  was  constructed  from  the  stand 
anl  linear  isoparametric  solid  element  with  a specified  bulk  modulus  and  zero  shear  modu- 
lus; and,  in  place  of  the  standard  2x2  t.aussian  integration  scheme,  one  point  at  the  ele- 
ment center  was  used. 

When  this  element  was  incorporated  into  the  code,  another  attempt  was  made  at  exe- 
cuting the  axisymmetrie  model;  this  proved  successful  and  verified  the  viability  of  the  ele- 
ment No  distortional  components  were  transferred  into  the  interior  elements,  and  the  ele- 
ment pressures  compared  favorably  with  a hydrostatic  distribution  (pressures  varied  from 
3 4 to  4 0 psi  throughout  the  interior  elements). 


N K I in  Discussion  ol  l*a pci  No.  751  1 f» I'lic*  development  ol  .1  detailed  finite  element  brain  model. " s,\l  Trans 
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I lie  applied  load  used  in  the  finite  difference  study  was  not  typical  ol  skull  impact 
loads  Originally,  fngin  applied  a step  pulse;  for  comparison,  Merchant  and  (irispino  chose 
the  same  load  magnitude  as  fngin.  However,  thc\  terminated  the  pulse  at  0 0614  ms,  as 
shown  m figure  4 13  which  was  lar  below  typical  load  pulse  durations  ( I he  initial  disturb 
ancc  reached  onh  two-thirds  of  the  sphere  diameter  by  the  t ime  the  pulse  was  terminated. ) 
Nevertheless,  this  load  was  used  in  the  present  comparison  study.  It  was  applied  uniformly 
over  a half  angle  of  7.5  degrees  (equivalent  to  0 5-4  square  inches  of  surface  area)  and  was 
directed  radially  inward  (see  f igure  4 12). 

Pressure  history  comparisons  are  shown  at  two  points  (sec  I igurc  4 12)  along  the  axis 
colinear  with  the  load  I he  comparisons  of  results  arc  shown  in  I igurcs  4 14  and  4-15.  1 he 
time  step  si/c  employed  was  3 n s and  is  approximately  twice  the  critical  time  step  si/c  re- 
quired tor  the  finite  difference  study  i inite  element  data  shown  has  been  sampled  ever\ 
three  time  steps. 

Surprisingly  good  overall  agreement  can  be  seen  lor  the  pressure  time  history  within  the 
brain  element  nearest  the  pole  (f  igure  4-14).  I he  primary  difference  appers  to  be  a time- 
shift  or  lag  between  the  two  responses  Also  noted  is  a difference  between  the  responses  for 
the  peak  compressive  atiel  peak  tensile  pressures.  I'hesc  differences  are  minor  considering 
the  dissimilar  manner  in  which  the  skull  bone  was  modeled  and  the  dissimilarities  in  the 
finite  element  and  finite  difference  techniques. 

Pressure  magnitude  decreases  considerably  for  points  away  from  the  pole  as  show  n for 
point  2 (f  igure  4 15).  At  the  same  time,  agreement  between  the  two  responses  deteriorates. 
I wo  explanations  are  offered  beyond  those  already  mentioned.  I irst  is  that  perhaps  the 
coarse  sampling  rate  of  the  finite  element  data  prevented  better  agreement.  Second,  and 
probably  more  significant,  is  the  effect  of  the  two  areas  of  ill-conditioned  finite  elements 
along  the  / axis  (see  f igure  4-12).  In  the  proximity  of  these  elements,  finite  element  pres- 
sure v alues  appear  to  be  amplified  with  what  must  be  regarded  as  numerical  noise.  I lie  fi- 
nite difference  mesh  was  better  conditioned  and  did  not  exhibit  these  errors.  It  is  therefore 
believed  that  by  increasing  the  sampling  rate  of  data  and  by  improv  ing  the  mesh  the  agree- 
ment would  be  improved. 

With  regard  to  computation  of  wave  speeds  in  the  spherical  model,  finite  element  data 
agrees  well  with  wave  propagation  speeds  computed  by  finite  difference  and  In  analytical 
methods  for  example,  taking  a diametrical  path  through  the  brain  material,  the  first  cl  is 
turbance  at  the  counterpole  is  expected  at  about  0. 1 1 ms.  I aking  a circumferential  path 
through  the  skull  bone  material,  the  first  disturbance  should  arrive  at  the  counterpolc  in 
about  0.09  ms  Results  ot  the  finite  element  model  show  that  both  the  brain  wave  and  skull 
wave  have  arrived  at  the  counterpole  at  0.12  ms. 

I o study  the  influence  of  time  step  si/e  on  computed  dynamic  response,  the  a.xisvm- 
mctric  model  described  above  was  employed  A more  realistic  load  was  applied  so  that  the 
model’s  performance  is  more  properly  evaluated.  Patrick,  et  al.55  measured  force  histories 


I VI  I'atriek.  II.  J Mert/,  Jr.,  and  < K Kroell.  '‘Cadaver.  knee,  chest,  and  head  impact  loads,"  in  I'roeeedmjjs  of  the 
I levcnth  Ntapp  (art  rash  ( onference,  ( )ct  10  11,1 967. 
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from  skull  collisions  against  padded  force  transducers,  f rom  this  data  a force  of  2,650 
pounds  with  a duration  of  10  ms  was  selected  as  a typical  head  injury  load.  An  equivalent 
pressure  (distributed  over  0.530  sq  in.)  was  computed  anil  applied  to  the  model.  I his  equiv- 
alent pressure-time  historv  is  shown  in  f igure  4-16. 

Results  of  this  studs  were  encouraging  for  two  reasons,  first  of  all  working  with  typical 
load  pulse  durations,  the  computed  finite  element  data  was  found  to  be  smooth  anti  not  ob- 
scured bv  jagged  response  data  as  was  seen  tor  the  short  duration  pulse.  Second,  accurate 
results  were  obtained  with  very  large  time  steps,  and  thus  analyses  of  typical  load  pulses  are 
expected  to  be  economical  and  suited  to  the  finite  element  method. 

Since  no  changes  in  the  axisymmetric  model  were  made,  the  critical  time  step  size  re- 
mained the  same,  about  0.0026  ms.  Computer  runs  were  made  for  each  of  the  following 
specified  time  step  sizes  0.005,  0.1 , 0.2,  and  0.4  ms.  Thus,  ratio  values  of  specified  to 
critical  time  step  size  of  about  2,  40,  80,  and  160  were  investigated. 

f igure  4 17  shows  the  compressiv  e stress  in  the  outer  table  bone  layer  directly  beneath 
the  site  ol  the  applied  load  No  noticeable  difference  in  the  stress  response  is  observed 
among  the  investigated  time  step  size  parameters.  Also  the  results  indicate  that  the  model 
satisfies  equilibrium,  f igure  4 18  shows  the  pressure  history  response  within  the  brain  at 
points  1 and  2 (see  f igure  4 12)  Again  no  noticeable  difference  is  observed.  These  results 
demonstrate  that  no  accuracy  loss  occurs  due  to  the  use  of  very  large  time  step  sizes. 

Although  in  this  study  the  main  intent  was  to  demonstrate  only  the  applicability  of 
finite  clement  modeling  to  head  injury  study,  some  data  pertinent  to  head  injury  was  col- 
lected and  is  shown  in  f igure  4 19.  Stress  profile  data  were  collected  from  the  computer 
run  associated  with  the  largest  specified  time  step  size,  0.4  ms.  The  buildup  and  decay  of 
compressive  and  tensile  pressures  w ithin  the  brain  at  the  pole  and  counterpole,  respectively, 
can  be  observed.  Note  that  there  is  no  indication  that  tensile  pressures  exist  at  the  poles 
directh  beneath  flic  applied  load.  I his  is  contrary  to  the  data  shown  previously  (Figure 
4 X1  and  is  ascribablc  to  the  relatively  long,  but  more  typical,  duration  of  the  load  pulse 
applied  :n  this  study 

l o verity  the  predicted  or  computed  rigid  body  motion  of  the  axisymmetric  model,  the 
theoretical  rigid  body  displacement  solution  was  first  generated  in  the  manner  of  particle 
dynamics.  Knowing  the  composite  mass  properties  of  the  model  (Table  4 1 ) and  the  forcing 
function  applied  to  the  model  (f  igure  4 15),  a closed  form  expression  for  the  acceleration 
can  be  written  and  then  integrated  twice. 

I he  results  of  the  comparison  are  shown  in  f igure  4-20.  It  can  be  seen  that  the  finite- 
element  data  is  in  good  agreement  with  the  theoretical  rigid  body  displacement  solution  at 
those  times  when  the  applied  load  is  nearly  terminated  (t>0.0()8  second).  Prior  to  that  time 
the  periodical  displacement  component  is  significant,  and  since  averaging  finite  element  data 
tends  to  nullify  this  component  the  agreement  can’t  be  expected  to  be  as  good.  Alternative- 
ly, neglecting  the  periodical  component  in  the  theoretical  solution  produces  good  agreement 
at  all  times. 

I he  foregoing  results  from  the  axisymmetric  model  enable  the  following  conclusions 
regarding  the  applicability  of  the  finite  element  modeling  technique  in  head  injury  studies. 

( 1 ) Acceptable  agreement  between  finite  element  and  finite  difference  results  for  an 
axisymmetric  head  injury  model  was  achieved  for  a short  duration  load  pulse. 
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(2)  Specified  step  sizes  associated  with  the  implicit  integration  scheme  of  about  2,  40. 
HO,  and  160  times  critical  time  step  size  exhibited  nearly  identical  results.  This  is  indicative 
of  the  efficient  utility  of  the  implicit  scheme  in  head  injury  studies  with  typical  applied  load 
durations. 

(3)  I he  time  history  of  rigid  body  translation  was  predicted  accurately  for  the  axisym- 
metrie  head  injury  model. 

(4)  Simulation  of  brain  material  was  achieved  successfully  be  degenerating  the  isopara- 
metric element's  spatial  integration  to  one-point  quadrature  at  the  element  center,  and  by 
specifying  the  shear  modulus  and  the  bulk  modulus  value  for  brain  material. 

(5)  l o the  extent  that  an  axisymmetric,  spherical  shell  is  a valid  head  injury  model,  the 
finite  element  method  is  shown  to  be  an  accurate  method  for  head  injur\  analysis. 

(6)  These  results  and  the  versatility  of  the  finite  element  method  encourages  its  appli- 
cation to  more  geometrically  complex  head  injury  models. 

Plane  Strain  Model 

Of  interest  in  this  model  study  was  the  influence  of  more  realistic  geometry  on  skull 
bending  and  on  distribution  of  intracranial  pressures.  Knowledge  acquired  of  these  para- 
meters could  contribute  to  what  has  been  learned  previously  concerning  head  injury  from 
axisymmetric  models  and  also  intimate  what  can  be  expected  in  the  eventual  extension  to 
three-dimensional  modeling. 

A plane  strain  model  was  constructed  to  simulate  the  geometry  of  a unit  slice  of  the 
midsagittal  plane  of  the  human  cranium.  Phis  model  is  shown  in  figure  4-2 1 . I he  geom- 
etry was  developed  from  measured  data  taken  from  the  plane  of  symmetry  of  a plastic  rep- 
lica of  the  human  skull.  A total  of  420  four-node  quadrilateral  elements  and  469  nodes  con- 
stitute the  discretization.  As  in  the  axisymmetric  model,  the  skull  bone  is  simulated  with 
three  layers  of  elements  through  the  thickness,  and  the  reduced  quadrature  elements  charac- 
terize the  encapsulated  brain  matter.  Loading  was  prescribed  as  a uniformly  varying  pres- 
sure with  a variation  in  the  form  of  a haversine  function.  I he  duration  of  loading  was  10  ms 
which  is  typical  of  durations  in  vehicle  accidents.  A time  step  size  of  0.4  ms  was  employed 
in  the  integration  ol  the  equations  of  motion.  Boundary  conditions  were  chosen  to  simu 
late  restraint  against  large  rigid  body  motions  at  the  base  of  the  cranium  near  the  neck  junc- 
ture. In  this  model  the  skull  was  treated  as  a closed  container. 

Previous  experimental  measurements  of  pressures  in  the  midsagittal  plane  have  indicated 
pressure  gradients  to  be  linearly  varying  from  compression  to  tension  in  an  anterior  to  pos- 
terior (A-P)  sense  Results  of  the  plane  strain  model  demonstrate  that  this  could  only  be 
true  for  pure  translation  of  the  skull  in  the  A-P  direction.  In  actuality,  for  the  same  blow, 
the  pressure  profile  will  be  determined  by  the  neck  constraints  and  the  resulting  relative 
"mixture”  between  translation  and  rotation  of  the  head.  In  figure  4 22  the  pressure  pro 
liles  are  quadratic  in  the  A-P  direction  with  considerably  more  tensile  pressure  than  would 
result  from  a linear  profile,  f igures  4 23  and  4-24  show  the  pressure  histories  in  the  brain 
and  the  stress  histories  in  the  outer  table  bone,  respectively. 

I wo  additional  computer  runs  were  made  with  the  two-dimensional  plane  strain  model 
ot  the  midsagittal  plane,  flic  objective  of  these  two  runs  was  to  evaluate  the  influence  of  a 
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viscoelastic  shear  modulus  characterization  of  brain  matter.  The  geometry  possessed  a 
“vented,”  instead  of  closed,  skull  configuration  by  providing  an  opening  at  the  foramen 
magnum.  Constraints  were  provided  against  translation  and  rotation  near  the  opening,  and 
the  dynamic  load  was  applied  to  the  forehead  region  just  as  for  the  previous  runs  with  this 
model.  I he  material  properties  for  the  brain  shown  in  I able  4-2  were  chosen  so  that  the 
shear  modulus  was  strongly  viscoelastic,  thus,  the  initial  shear  modulus  value  was  1 5,500 
psi  and  at  the  end  of  10  ms  the  shear  modulus  decayed  to  2,242  psi.  While  the  shear  re- 
sponse was  made  strongly  dependent  on  time,  the  bulk  modulus  was  considered  as  elastic. 

I o establish  a basis  for  comparison,  a completely  elastic  characterization  of  brain  matter 
was  employed  in  the  second  -.omputer  run.  figure  4-25  shows  the  stress  histories  for  both 
computer  runs  Initially,  no  difference  in  the  response  was  noted,  but  eventually  the  inher- 
ent damping  properties  of  the  viscoelastic  material  manifest  themselves  by  dissipating  the 
energy  at  a continuously  increasing  rate  with  time.  At  7.5  ms,  for  example,  approximately 
one-half  the  elastic  stress  response  has  been  damped  by  the  viscoelastic  characterization, 
figure  4-26  shows  the  bending  stress  history  in  the  outer  fibers  of  the  skull  bone  at  the  site 
of  the  blow.  Since  the  brain  is,  in  effect,  softer  in  the  viscoelastic  case,  the  elastic  skull  bone 
must  deflect  more  and  develop  higher  stresses  than  the  elastic  brain  case. 

Obviously  a v iscoelastic  characterization  has  a great  influence.  I he  difficulty  is  in  the 
determination  of  a more  accurate  and  representative  viscoelastic  shear  modulus  formulation. 

I he  one  used  in  this  analysis,  shown  in  I able  4-2,  was  arbitrarily  based  on  a Poisson’s  ratio 
value  of  0.475;  but  from  this  information  an  idea  can  be  gained  as  to  the  sensitivity  of  dy- 
namic skull  response  to  viscoelastic  brain  matter  1 his  information  will  assist  in  the  charac- 
terization of  brain  material  properties  in  subsequent  three-dimensional  model  analvses. 

fhe  computed  bending  shape  of  the  skull  resulting  from  the  loading  and  boundary  con- 
ditions described  above  is  presented  in  f igure  4 27.  Since  this  is  a linear  analysis,  only  the 
shape  (as  opposed  to  magnitude)  is  of  interest,  fhe  largest  bending  stresses  occurred  as 
tensile  stresses  on  the  inside  cranial  surface  adjacent  to  the  location  of  applied  pressure. 

I ensile  stresses  are  known  to  be  conducive  to  linear  fracture  in  cranial  bone.  I hey  arc  indi- 
cated on  either  side  ol  the  skull  surface  where  they  occur.  Note  the  location  of  inflection 
points  which  indicate  a sign  reversal  for  bending  stresses  in  the  cranial  bone.  Thus,  while 
the  model  suggests  that  the  area  most  conducive  to  linear  skull  fracture  is  on  the  inside  sur- 
face near  the  loading,  an  alternative  location  would  exist  on  the  ouside  surface  near  the  top 
of  the  skull. 

All  bending  stresses  in  the  skull  were  observed  to  rise  and  subside  synchronously  with 
the  havcrsinc  loading  function  Displacement  components  due  to  bending,  therefore,  would 
also  follow  the  same  time  historv  However,  the  resultant  movement  lagged  behind  the 
loading  history.  In  this  instance  a verv  slight  counterclockw  ise  rotation  of  the  cranium  was 
observed  as  the  applied  pressure  displaced  it  “up  and  back  over”  the  restraining  support. 

I urthermore,  this  motion  occurred  while  the  applied  pressure  was  subsiding  and  was  there- 
fore, considered  to  be  an  inertial  effect. 

I his  behav  ior  is  further  evidenced  by  the  intracranial  pressure  data  shown  in  f igure 
4 28  Pressure  contours  continued  to  form  well  after  the  applied  load  has  peaked  and  begun 
to  subside  (t  • 5 ms).  Contour  values  normalized  on  the  applied  pressure  are  listed  in  fable 
4 3 Brain  damage  is  often  observed  clinically  to  occur  on  the  opposite  side  of  the  skull 


I able  4 2 brain  Material  Properties  for  Plane  Strain  Model 


< haracterization 

K ( psi) 

< i (psi) 

1 last  ie 

305,000 

15,500 

Viscoelastic 

305,000 

(l(t)a 

1 luid 

305,000 

0 

1 

aCi(t)  5,800  e‘l32t  + 9,700c  264t 

(.MM  15,500 

i'(0)  0 475 

(i(0.01 ) 2,242 

/'(().()  1 ) 0.4988 


I able  4 3 Normalized  Pressure  < Contour  Values 


Contour  Number 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 


Pressure  Ratio3, 

P7Pa 

0.350 
-0.265 
0.180 
-0.096 
0.01  1 
+0.073 
+0  158 
+0.243 
+0  328 
+0  413 


‘‘Minus  sign  indicates  compression  in  brain,  positive 
sign  indicates  tension  in  brain 


from  site  ot  the  blow  and  is  termed  contrecoup  damage  I vidence  indicates  that  lesions  oc 
curring  in  this  manner  on  the  surface  of  the  brain  are  attributable  to  negative  f tensile  > pres 
sures  which  rupture  surface  capillary  blood  vessels.  Determining  the  location  of  negative- 
pressures  is,  therefore,  an  important  function  of  a head  injury  model.  I'he  contours  show 
that  the  location  of  the  maximum  negative  pressure  occurred  near  the  back  and  top  of  the 
brain  1 his  result  is  related  to  the  direction  of  inertia  forces  which  are,  in  turn,  governed  b\ 
tlte  particular  combination  of  loading,  geometry,  and  restraints  employed  \ll  three  of 
these  parameters  must  be  simulated  accurately  in  head  injury  models  because,  as  seen  in  this 
example,  they  have  an  important  influence  on  predicting  both  skull  bending  modes  and 
intracranial  pressure  distributions. 
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I'he  two  dimensional  models  discussed  in  the  preceding  chapter  were  intended  to  he 
preparaton  tor  a three-dimensional  anal)  sis.  Since  the  finite  element  method  is  a viable 
three-dimensional  anal)  sis  technique,  there  seemed  little  reason  to  dwell  on  the  two- 
dimensional  studies  after  thev  had  served  their  preliminary  purpose,  and  to  proceed  toward 
the  ultimate  goal  of  the  project  a three-dimensional  model.  It  is  felt,  however,  th.n  these 
simple  models  can  be  exploited  further  that  there  is  more  to  learn  troni  them  and  should 
not  he  discarded. 

\ primarv  consideration  in  the  develop  tent  of  the  HIM  code  was  attainment  ol  accur- 
ate geometrical  simulation  of  the  cranial  anatomy.  I he  term  “recognizable  geometry”  is 
employed  to  distinguish  this  feature  from  previous  head  injury  model  studies  employing 
rotationalh  s\  mmctric  geometry. 

I'he  effort  to  develop  a recognizable  three-dimensional  discretization  of  the  skull  paral- 
leled the  preliminary  model  study  effort  as  pointed  out  in  the  OKCiANI/A  I K ) N HI  AN 
section  2 As  expected,  it  was  a very  involved  procedure,  taking  approximately  1 t ear  from 
start  to  finish.  I his  chapter  describes  the  procedure  and  the  considerations  made  during 
that  time  I his  description  should  familiarize  potential  users  of  the  HIM  code  with  various 
parts  of  the  total  discretization,  the  various  orders  of  discretization,  and  other  options  bear- 
ing on  effective  use  of  the  model. 

( omponents  of  the  Discretization 

I'he  first  consideration  in  discretizing  the  skull  was  to  determine  which  anatomical  struc- 
tures should  be  included  Besides  the  cranial  bone  and  brain  structures  (obvious  compo- 
nents of  a head  injury  model),  other  structures  could  have  required  inclusion,  \mong  these 
w ere  the  neck,  facial  bones,  subarachnoid  space,  and  the  major  intracranial  membranes  such 
as  the  falx  and  tentorium.  I he  inclusion  of  these  structures  and  other  more  intricate  intra- 
cranial structures  had  to  be  weighed  with  the  added  complexity1  and  cost  associated  with 
meaningful  discretizations  of  each.  Their  relative  importance  had  to  be  assessed  In  the  ex- 
tent to  which  each  participates  in  common  head  injury  mechanisms  l him. itch  , the  decid 
ing  factor,  however,  had  to  be  whether  (lor  feasibility  or  economic  reasons)  certain  of  these 
structures  could  be  included  in  a finite  clement  model  I'he  result  of  having  omitted  some 
of  these  structures  . ould  have  adversely  affected  tin  scope  of  the  model's  applicability  ; to 
obtain  reliable  information,  the  effort  nevertheless  had  to  proceed. 

Skull  Bone  Struct  ure  \ detailed  study  of  the  anatomy  of  skull  sutures  was  made  with 
the  intent  of  incorporating  their  influence  on  the  structural  behavior  of  the  skull.  I ight 
bones  constitute  the  human  cranium  frontal,  occipital,  two  temporal,  two  parietal  and  two 
sphenoidal  as  shown  in  I igure  5 1 I he  sutures  shown  in  f igure  5-2  are  the  natural  inter 
faces  found  between  these  bones  and  may  form  a failure  mechanism  under  certain  loading 
conditions  and  for  certain  age  groups  \t  birth,  unossified  membranous  intervals  (fontanels) 
exist  between  the  skull  bones  but  will  close  within  two  month  after  birth  except  the  anter- 
ior which  closes  at  about  21  .>  years  of  age.  Apparently  this  process  merely  closes  the  gap 
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between  bones  but  Joes  not  solidify  them  together.  Suture  closure  begins  at  about  age  22 
\ ears  when  ossification  converts  the  soft  hyaline  cartilage  coating  at  the  joints  into  bone, 
first  the  sagittal  and  sphenoidal  fontanels  ossify,  then  at  age  24  years  the  coronal  suture 
ossifies,  and  at  age  26  the  lambdoidal  and  occipitomastoid  sutures  ossify.  Between  the  age 
of  26  and  30,  ossification  is  very  rapid  compared  to  thereafter  when  it  slows  down.  A final 
burst  of  ossification  occurs  in  old  age  Apparently  this  ossification  is  independent  of  sex. 

A trial  discretization  of  the  skull  surface  which  included  “suture  elements”  was  made, 
anticipating  the  specification  of  mechanical  properties  for  these  elements  (see  f igure  5-3  ). 
Difficulty  had  been  envisioned  when  the  entire  cranium,  including  the  contents,  was  to  be 
discretized  in  a manner  consistent  with  the  suture  discretization ; the  two  seem  to  be  con- 
flicting requirements.  1 hus,  the  idea  of  including  suture  elements  was  discarded.  At  the 
same  time  the  scope  of  the  head  injury  model  might  be  more  applicable  to  adults  than  to 
younger  individuals. 

Studies  with  a previous  finite  element  model  of  the  cranial  structure18  concluded  that 
the  structural  significance  of  the  sandwiched  or  layered  construction  existing  naturally  in 
the  skull  bone  is  important  and  should  be  included  somehow  in  the  simulation  Though 
this  work  made  use  of  composite  shell-type  elements  and  the  present  effort  will  make  use  of 
solid  elements,  it  is  a simple  matter  to  effect  this  simulation.  A decision  to  include  a layered 
skull  discretization  in  the  HIM  code  was  primarily  based  upon  the  above  cited  recommenda 
tion.  This  recommendation  was  based  upon  consideration  for  proper  skull  strength  Simula 
tion.  Yet  there  was  still  another  advantage  in  the  present  investigation  for  doing  this. 

Actual  skull  bone  is  comprised  of  three  distinct  layers  of  bone  through  its  thickness.  The 
innermost  and  outermost  layers  are  a dense  bone  material  and  are  referred  to  as  the  inner 
and  outer  table  bone.  The  middle  layer  is  by  contrast  a very  porous  bone  material  and  in 
general,  is  thicker  than  the  other  two  layers;  this  layer  is  referred  to  as  the  diploe  layer  By 
simulating  this  structure  with  three  layers  of  solid  elements,  one  not  only  captures  the  na- 
tural variation  of  skull  bone  material  stiffness  through  the  thickness  (transverse  anistropy), 
but  also  presumes  to  obtain  the  actual  stress-strain  behavior  of  the  skull  bone.  I his  is  im- 
portant because  such  information  can  be  correlated  directly  with  skull  bone  fracture. 
Whereas  a simulation  with  a single  layer  of  shell  elements  can  adequately  represent  the  over- 
all stiffness  of  the  skull  bone,  it  cannot  yield  useful  stress-strain  information  because  the 
assigned  material  properties  must  of  necessity  be  averaged  or  smoothed  values. 

I here  are  several  reasons  why  the  cranial  structure  should  be  simulated  carefully  despite 
the  fact  that,  generally,  brain  damage  can  occur  from  impact  in  the  absence  of  skull  frac 
ture.  A linear  head  injurs  model  must  be  capable  of  predicting  the  onset  of  skull  fracture, 
if  not  the  fracture  process  itself  Inclusion  of  the  skull  bone  renders  the  model  more  usable 
in  the  management  of  brain  trauma  When  brain  damage  accompanies  a severe  fracture,  the 
prognosis  can  be  considerably  complicated  by  intracranial  infection. 

The  elastic  and  inertial  forces  causing  brain  tissue  alteration  are  partially  the  result  of 
the  particular  way  in  which  the  skull  deforms  and  absorbs  impact  energy  in  each  instance 
I hus,  notwithstanding  that  skull  fracture  is  in  itself  an  important  component  of  injury,  cor 
rect  simulation  of  the  cranial  bone  structure  is  also  necessary  in  order  to  predict  the  input 
forces  to  the  brain 
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Linear  skull  fractures  are  the  most  commonly  observed  types  of  skull  fracture.  Linear 
fractures  tend  to  occur  radially  front  a point  of  maximum  deflection  but  this  pattern  is 
somewhat  modified  by  irregular  areas  of  bony  strength.  Two  observations  are  important  in 
regard  to  linear  fracture.  First,  when  the  blows  cause  very  localized  bending  the  inner  table 
is  more  severely  fractured  than  the  outer  table.  Second,  the  location  of  the  blow  and  the 
path  followed  by  the  propagating  fracture  are  related.  l:or  frontal  blows  the  paths  tend  to 
radiate  and  very  often  fracture  the  thin  bone  of  the  orbital  roofs.  For  vertical  blows  the 
paths  classically  radiate  outward  over  the  top  of  the  skull.  For  occipital  blows,  the  paths 
run  toward  the  base  and  the  foramen  magnum,  f inally,  for  parietal  blows  the  paths  radiate 
toward  the  temporal  bone.  Depressed  skull  fractures  are  usually  caused  by  more  blunt  im- 
pacts and  are  characterized  by  depression  or  collapse  of  the  cranial  bone’s  sandwiched  struc- 
ture. Large  flat  areas  of  impact  can  cause  extensive  collapse  of  the  bone  and  greatly  modify 
the  subsequent  energy  transmitted  to  the  brain. 

The  skull  bone  structure  can  be  thought  of  as  comprised  of  cranial  bones  and  facial 
bones.  I he  cranial  bone  complex  contains  and  protects  the  brain  and  is  the  more  important 
of  the  two  types.  T he  skull  bone  discretization  should  reflect  this  importance,  but  the  in- 
clusion of  the  facial  bones  is  necessary  because  their  inertial  properties  substantially  influ- 
ence the  resulting  motion  (acceleration)  of  the  skull.  As  it  happens  in  finite  element  tech- 
nology, a reasonable  simulation  of  a structure’s  inertial  properties  can  be  effected  with  a 
surprisingly  coarse  discretization.  (The  same  is  not  true  of  a structure’s  elastic  properties.) 

In  this  case  all  that  is  required  is  a reasonable  simulation  of  the  inertial  properties  while 
no  quantitative  stress-strain  information  is  sought  in  the  facial  bone  structure.  T herefore, 
a coarse  facial  bone  discretization  has  been  included  in  the  HIM  code  geometry. 

Another  important  consideration  is  the  neck  structure  and  whether  or  not  it  should  be 
included  as  part  of  the  head  injury  model.  What  is  of  primary  concern  is  that  many  so- 
called  head  injuries  are  actually  neck  injuries  or  neck -related  injuries;  for  example,  damage 
to  the  spinal  cord  and  brainstem  from  either  the  hyperextension  or  hyperflexion  modes  of 
injury  is  very  common.  A complete,  all-inclusive  head  injury  model  should  be  able  to  simu- 
late them,  f urthermore,  damage  to  the  brain  stem  even  though  caused  by  hyperextension 
or  hyperflc.vion  is  often  intracranial  damage  because  the  brainstem  is  largely  situated  within 
the  cranium. 

From  a practical  modeling  standpoint  however,  there  is  a jurisdictional  dilemma  arising 
from  the  discontinuous  nature  of  the  head/neek  skeletal  structure  and  the  continuous 
nature  of  the  central  nervous  system.  I he  anatomical  structure  of  the  vertebral  column  sug 
gests  strongly  that  the  neck  be  modeled  continuously  with  a thoracic  model56  instead  of 
with  a skull  model.  Yet,  because  the  brainstem  and  spinal  cord  are  continuous  at  the  head/ 
neck  interface,  deleting  the  neck  from  the  skull  model  is  undesirable  because  the  deletion 
limits  the  scope  of  the  head  injury  model  to  injury  prediction  exclusive  of  neck -related 
trauma  At  this  time  practical  economics  and  modeling  feasibility  prohibit  the  ideal  con- 
sideration ot  a combined  head/neck  model. 
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Without  a neck  structure  discretization  interacting  at  the  base  ol  th,  skull,  the  forcing 
function  of  the  neck  on  the  skull  base  must  be  treated  as  necessary  input  information  to  the 
skull  model.  1 his  force  is  extremely  difficult  to  estimate.  I his  basic  difficulty  is  not  re- 
moved with  the  inclusion  of  a good  neck  discretization  because  the  same  problem  would 
exist  with  the  specification  of  an  input  forcing  function  at  the  base  of  the  neck. 

brain  and  Intracranial  Structure.  Simulation  of  the  proper  support  conditions  for  the 
brain  is  ot  . indamental  importance  since  the  brain  is  not  rigidly  supported  within  its  con- 
iines (see  figure  5-4  for  reference  to  brain  anatomy).  It  is  tethered  at  its  base  by  the  exit 
ot  cranial  nerves  and  the  brainstem  from  the  cranium.  The  falx  and  tentorium  membranes 
lend  major  support.  1 hese  membranes  are  not  shown,  but  the  falx  lies  in  the  midsagittal 
plane  and  extends  downward  between  the  two  cerebral  hemispheres.  The  tentorium  forms 
a partition  between  the  cerebrum  and  cerebellum.  1 he  semisolid  medium  of  the  brain  is 
surrounded  bv  cerebral  spinal  fluid  (CST)  which  is  displaced  when  the  brain  shifts  from  one 
tethered  position  to  another. 

Intracranial  volume  shifts  induced  by  impact  give  rise  to  two  basic  types  of  injury. 

( oup  brain  injuries  occur  if  at  the  point  of  impact  the  energy  is  sufficient  to  damage  the 
brain  tissue.  I he  extent  of  damage  w ill  be  proportional  to  the  applied  force  magnitude. 

( ontrccoup  brain  injurv,  occurring  opposite  the  impact,  will  be  dependent  upon  the  direc- 
tion of  the  applied  force  vector  as  well  as  the  magnitude.  This  direction  will  determine 
w here,  along  the  skull/brain  interface,  the  brain  tends  to  separate  greatest  from  the  skull, 
flic  seriousness  of  the  injury  will  depend  on  whether  or  not  bony  irregularities  exist  or 
whether  existing  nerve  roots  exist  at  that  particular  section  of  cranial  wall. 

Three  kinds  of  brain  tissue  alteration  or  failure  can  be  defined  as  a result  of  impact. 

I hcv  are  concussion,  contusion  (bruise),  and  compression  of  the  brain.  The  first  is  not  ex- 
plained bv  clinical  or  pathological  findings  as  are  the  latter  two.  Concussion  is  attended  bv 
loss  of  consciousness,  no  matter  how  brief,  and  is  caused  either  by  damage  to  the  cortex  or 
to  the  arousal  mechanism  in  the  brainstem.  Contusions  of  the  brain  are  generally  associated 
w ith  the  coup  anil  contrecoup  iipury  mechanisms  and  arc  pathologically  observable,  as  is 
compression  of  the  brain  that  results  from  internal  bleeding  (subdural  hematomas).  Ana- 
tomically the  lesions  of  concussion,  contusion,  and  compression  are  distinctly  different, 
because  of  this  and  because  they  were  not  simply  related,  each  will  most  probable  require 
a separate  failure  theory.  There  is  no  reason  to  believe  that  such  theories  cannot  be  based 
upon  the  fundamental  mechanical  phenomenon  ot  stress  or  the  geometrical  pheonomenon 
of  strain  If  not  stress  and  strain,  then  one  might  ask,  what  other  factors  arc  there?  The 
quantities  of  velocity,  acceleration,  and  force  must  all  ultimately  manifest  themselves  in 
some  form  of  stress  and  strain  patterns  within  the  brain  tissue.  It  seems  as  though  these 
patterns  offer  the  best  and  most  direct  route  to  predicting  brain  damage. 

As  can  be  seen  in  f igure  5-4,  the  human  brain  is  structurally  very  complex.  It  becomes 
obvious  that  some  degree  of  idealization  is  necessary  when  discretizing  this  structure.  The 
analysis  tools  with  which  we  must  proceed  (finite  element  technology)  do  not  easily  allow 
separation  (displacement  discontinuities)  of  the  various  intracranial  components  depicted, 
even  though  they  arc  separate  components.  II  this  capability  were  available,  the  finite  clc 
ment  method  could  allow  an  analv  sis  which  would  include  these  components  as  separate 
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entities  in  the  bruin.  It  is  not  presently  available,  therefore,  the  major  idealization  in  model- 
ing till'  brain  is  the  consideration  of  a homogeneous  intracranial  structure  possessing  a geom- 
etry roughly  equivalent  to  that  of  the  intracranial  v ault. 

In  any  dynamic  finite  element  stress  analysis  the  two  internal  forces  sought  implicitly  are 
the  inertial  and  elastic  lorces.  The  inertial  force  is  merely  the  product  of  the  mass  and  accelera- 
tion of  a typical  node  point  in  the  finite  element  mesh.  II  the  individual  mass  densities  of  the 
intracranial  substructures  arc  fairlv  uniform  then  the  field  ol  inertial  forces  will  be  homogeneous 
and  everywhere  directlv  proportional  to  the  acceleration  field.  I herefore.  correct  simulation  ol 
the  inertial  force  field  is  not  significantly  affected  by  the  assumption  of  a homogeneous  intra- 
cranial structure  since  intracranial  materials  do  have  mass  densities  that  are  fairly  uniform. 

l'h is  is  not  the  case  with  internal  elastic  forces.  I hese  forces  are  dependent  upon  the 
correct  simulation  of  the  conditions  of  displacement  continuity  between  the  intracranial 
substructures.  I hese  conditions  are  not  always  continuous  and  are  difficult  to  define.  Of 
course,  the  elastic  forces  are  also  dependent  upon  the  mechanical  properties  of  these  mater- 
ials but  this  is  not  the  difficulty,  liven  though  these  mechanical  properties  arc  known  or 
can  be  estimated  fairlv  accurately  and  can  be  individually  prescribed,  the  elastic  forces  still 
cannot  be  accurately  simulated  unless  the  displacement  continuity  (or  discontinuity)  con- 
ditions can  first  be  correctly  prescribed.  Since  this  cannot  be  done  as  a practical  matter 
w ith  present  state-of-the-art  capability,  it  appears  inadvisable  to  think  in  terms  of  more  de- 
tail in  the  construction  of  the  brain  discretization.  In  this  connection,  it  should  be  men 
tinned  that  mesh  uniformity  has  import  beyond  that  of  esthetics.  Where  possible,  the  mesh 
discretization  should  be  constructed  unif  ormly  for  the  sake  of  computational  accuracy. 

I his  is  especially  true  for  certain  finite  element  types  such  as  those  based  upon  incompatible 
finite  element  formulations.  Naturally,  more  detail  in  a discretization  puts  a needless  strain 
on  the  uniformity  requirement  and  will  adversely  influence  numerical  accuracy  in  those 
cases  where  the  additional  detail  is  ineffectual  or  unwarranted. 

I here  are  some  exceptions  to  the  idealization  approach  which  utilizes  a completely 
homogeneous  intracranial  structure  and  uniform  discretization.  These  will  be  discussed 
separately  in  sections  following. 

Since  the  brain  is  not  attached  continuously  along  the  internal  cranial  wall  but  is  sup 
ported  and  tethered  at  their  interface,  a method  must  be  devised  to  approximate  this  sus 
pension.  The  single  layer  of  elements  which  simulate  the  subarachnoid  space  was  prov  ided 
in  the  basic  HIM  code  configuration  for  this  eventuality.  In  this  wav  provision  is  made  for 
the  possibility  of  relative  motion*  between  the  skull  and  brain  at  their  interface.  The  ana 
tomy  suggests  that  this  motion  might  occur  because  the  brain  is  not  continuous  with  the 
skull  Indeed  the  interface  appears  to  be  lubricated  with  cerebral  spinal  fluid  Also  neuro- 
surgeons speak  a great  deal  about  intracranial  volume  shifts  and  the  brain  slosh” 
mechanism3,  where  the  primate  brain  has  been  observed  through  Incite  calvariums  to  dis- 
place relative  to  the  internal  cranial  wall,  but  just  providing  a layer  of  elements  does  not  in 


*ln  this  context , relative  motion  means  relative  motion  in  both  a normal  and  a tangential  direction  at  a point  on  the  skull 
f>raui  surface*. 

5"l<  II  I’udenz  and  ( II  Sheldon  I he  lueile  ealvarium  a method  for  direct  observation  of  the  brain.  Journal  ot 
Neurosurgery,  vol  3.  I'Mft,  pp  4H7  sos 
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itsclt  simulate  the  subarachnoid  space.  An  appropriate  set  of  material  properties  must  be 
assigned  to  those  elements,  depending  on  what  the  “physics”  of  the  situation  is  believed 
to  be.  1 he  particular  set  of  mechanical  properties  recommended  for  use  in  the  HIM  code 
is  discussed  in  Chapters  6 and  7,  but  it  should  be  said  here  that  this  specification,  like  any 
other  material  property  specification,  is  optional  to  the  HIM  code.  The  user  may  prescribe 
any  set  of  material  property  values  for  any  region  of  the  discretization. 

I he  falx  and  tentorium  membranes  may  well  present  another  exception  to  the  concept 
of  a completely  homogeneous  idealization  of  the  intracranial  structure.  These  membranes 
separate  the  intracranial  vault  into  separate  compartments  and  can  provide  maj  r support 
which  would  resist  intracranial  volume  shifts.  This  is  an  important  structural  influence  and 
therefore  simulation  of  these  membranes  was  given  careful  consideration  in  the  HIM  code 
development,  besides  the  anatomical  evidence  supporting  a decision  to  include  them,  it  ini- 
tially appeared  a structual  simulation  could  be  easily  effected  with  basic  two-dimensional, 
membrane  finite  elements.  There  was,  however,  a dissenting  note  represented  by  the  experi- 
mental work  reported  by  Roberts,  et  al.58.  They  concluded  that  the  existence  of  the  mem- 
branes made  little  difference  to  the  intracranial  pressures  which  they  measured.  Neverthe- 
less, it  was  decided  initially  to  incorporate  the  effect  of  the  falx  and  tentorium  membranes 
into  the  HIM  code. 

Many  obstacles  appeared  in  the  initial  attempts  to  discretize  the  intracranial  contents. 
During  the  development  of  the  discretizing  process,  some  simplifications  were  necessarily 
made  to  the  model  to  achieve  a working  discrctizaiton.  One  of  these  simplifications  was 
deletion  of  the  requirement  for  the  membranes  and  relegation  of  their  inclusions  as  a rec- 
ommended feature  for  future  work.  I hus,  the  present  HIM  code  does  not  possess  the  influ- 
ence of  the  intracranial  membranes.  It  is  believed,  however,  that  the  mesh  discretizing  tech- 
nique used  can  eventually  be  modified  to  include  the  membranes. 

One  other  problem  relating  to  the  membranes  will  still  have  to  be  overcome.  It  is  not 
clear  how  the  two-dimensional  membrane  elements  can  be  kinetieallv  isolated  from  adjacent 
brain  elements  so  as  to  adequately  simulate  their  interaction.  As  these  membranes  exist  in 
the  skull,  they  are  not  affixed  continuously  with  adjacent  cerebral  material.  Tor  example, 
motion  of  the  cerebrum  in  an  anterior  posterior  direction  is  not  resisted  by  the  in-plane 
stiffness  of  the  falx  membrane.  However,  if  in  the  simulation,  the  membrane  elements  are 
not  correctly  isolated  (displacement  compatibility  relinquished),  the  adjacent  cerebral  ele- 
ments would  be  stiffened  artificially  and  incorrectly  in  the  anterior  posterior  direction. 

I herefore,  an  adequate  simulation  of  the  falx  and  tentorium  is  not  necessarily  achieved  by 
merely  including  two-dimensional  finite  elements  in  the  discretization  unless  provision  is 
also  made  for  their  proper  interaction  with  adjacent  three-dimensional,  cerebral  elements. 

It  is  believed  that  there  is  yet  a third  exception  to  the  homogeneous  intracranial  con- 
cept and  that  it  involves  the  brainstem  simulation.  Here  again  the  HIM  code  does  not  ac- 
count for  the  rather  fibrous  nature  of  the  brainstem’s  constitutive  makeup,  but  rather  as- 
sumes it  to  be  a homogeneous  material  identical  to  the  surrounding  intracranial  material. 


SH 

\ 1 Kohcrts,  V K Hodgson,  am I I M I homas  ' l-lind  pressure  gradients  caused  l»y  impact  to  thi  human  skull," 
Biomechanics  Monograph,  ASMI  . 1967 
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I lie  reasons  for  this  are  the  same  as  those  offered  in  the  discussion  of  the  membranes.  I irst, 
the  constraints  imposed  by  the  brainstem’s  geometrical  complexity  prevented  the  realization 
of  a viable  intracranial  discretization  and,  second,  the  inability  to  simulate  the  brainstem’s 
kinematic  interaction  with  surrounding  material  makes  any  attempt  to  do  so  inadvisable 
during  the  initial  development  of  a head  injury  model,  because  the  brainstem  is  generally 
thought  to  participate  widely  in  closed  brain  injury,  however,  it  does  deserve  additional  atten- 
tion  in  the  modeling  process. 

Methods  ot  Discretization 

At  the  outset  of  the  project,  geometrical  detail  of  the  model  was  assigned  major  impor- 
tance Automatic  mesh  generation  schemes,  which  are  frequently  associated  with  ordinary 
finite  element  problems,  were  initially  ruled  out  because  they  were  thought  to  be  efficient 
only  with  simple  or  analytically  describable  geometry.  A more  arduous  manual  approach, 
but  one  that  was  expected  to  be  more  appropriate  for  discretizing  complex  geometry,  was 
tirst  embarked  upon.  I his  decision  resulted  in  a delay  in  the  realization  of  a workable  finite 
element  discretization  of  the  skull-brain  system.  Though  eventually  an  automatic  mesh 
generation  scheme  was  employed  successfully,  it  is  believed  useful  to  first  discuss  the  origi- 
nal discretization  attempts. 

Manual  Mesh  (, dictation.  A life-sized  plastic  skull  replica  was  immersed  in  an  epoxy 
bath  for  subsequent  mechanical  discretization.  Once  the  assembly  was  completely  cured, 
it  w as  removed  from  the  mold  and  then  rectangularly  shaped  in  a flycutting  operation  anti 
thereafter  polished.  The  resulting  processed  cast  is  shown  in  figure  5-5.  Next  the  assembly 
was  sliced  into  thin  sagittal  sections,  each  ot  which  would  be  discretized  by  superimposing 
a grid  over  the  slices.  Measurements  relative  to  the  reference  coordinate  system  show  n would 
be  made  of  the  discretization  and  eventually  coded  tor  input  to  the  finite  element  model  ot 
head  I he  reference  coordinate  system  is  an  established  system  for  cranial  anatomy.  The 
\ and  v -axes  lie  in  the  I rankfort  plane.  I his  plane  contains  the  auditory  meatuses  and  the 
orbital  notches.  I he  v axis  is  co-lincar  with  the  line  connecting  the  auditory  meatuses,  and 
the  z-axis  is  formed  by  the  right  hand  rule. 

The  epoxv  cast  w as  sliced  four  times  in  planes  parallel  with  the  midsagittal  plane,  crcat 
ing  five  epoxv  slabs  The  slabs  are  shown  assembled  in  f igure  5-6.  A finite  element  mesh 
was  designed  and  superimposed  on  each  slab  f.ach  superimposed  mesh  was  composed  ex- 
clusively of  quadrilateral  subregions.  In  this  way  the  three-dimensional  elements  are  formed 
as  irregular,  six  sided  bricks,  and  the  discretization  is  devoid  of  tetrahedrons  or  triangularly 
shaped  elements.  In  addition  to  merely  discretizing  the  cerebral  volume,  the  finite  element 
mesh  has  been  designed  to  include  the  geometrical  detail  attending  the  falx  of  cerebrum 
membrane  and  the  tentorium  of  cerebellum  membrane. 

During  this  development  period  it  was  also  necessary  to  determine  the  number  ot  ele- 
ments and  nodes  required  for  an  adequate  structural  idealization  while  keeping  the  model’s 
size  economically  tractable.  Roughly  400  elements  constituted  the  initial  discretization 
which  included  about  375  brick  type  elements  ami  about  25  membrane  type  elements,  l or 
simulated  dynamic  load  applications  in  the  midsagittal  plane,  symmetry  was  taken  advantage 
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of  and  the  discretization  was  made  about  twice  as  fine.  However,  for  unsymmctrical  load 
application  the  coarseness  shown  is  approximately  what  must  be  used  because  400  addi- 
tional elements  will  have  to  be  generated  for  the  opposite  half  of  the  skull.  It  was  estimated 
that  800  to  1,000  elements  are  nearly  the  computer’s  storage  capacity  and  will  adequately 
represent  the  geometry  of  the  skull. 

1 he  difficulty  with  the  epoxy  cast  procedure  arises  from  the  parallel  plane  concept 
upon  which  it  is  based.  With  the  restriction  that  no  triangularly  shaped  elements  are  to  be 
allowed,  an  increasingly  finer  discretization  occurs  in  planes  farther  from  the  midsagittal 
plane.  This  presents  an  unacceptable  mesh  size  gradation.  It  would  be  better  from  the  stand- 
point of  mesh  gradation  if  the  fineness  of  discretization  would  increase  with  radial  distance 
inward  from  the  skull  surface  towards  the  center  of  the  brain,  because  with  this  latter  con- 
cept the  primary'  generators  of  the  mesh  must  be  radial  lines,  as  contrasted  with  parallel 
sections,  the  epoxy  east  procedure  would  have  no  particular  advantage.  Instead,  it  was 
thought  that  a tractable  approach  could  be  effected  with  pliable,  fine  wire  lengths  placed 
radially  in  the  plastic  skull  to  act  as  main  generators  of  the  mesh.  Two  difficulties  were 
foreseen  with  this  approach,  first,  a few  tetrahedron  elements  must  inevitably  be  included 
due  to  the  convergence  of  the  radial  generators.  Mixing  these  elements  with  six-sided  brick 
elements  was  to  be  avoided  for  computational  accuracy.  Second,  node  point  coordinate 
measurements  would  not  be  accomplished  as  easily  as  with  the  plane  slices  of  epoxy. 

Although  some  effort  was  expended  on  the  radial  wire  approach  and  the  original  epoxy 
cast  procedure  did  culminate  in  a usable  three-dimensional  discretization,  neither  method 
proved  to  be  satisfactory.  These  point-by-point  discretization  methods  were  too  arduous 
and  did  not  produce  sufficiently  uniform  discretizations  of  the  skull-brain  system.  As  a 
result  they  were  finally  abandoned  in  favor  of  an  automatic  mesh  generation  approach. 

Automatic  Mesh  ( feneration.  This  approach  to  the  discretization  of  the  skull-brain  sys- 
tem is  primarily,  but  not  completely,  automatic  in  its  computation  of  nodal  point  coordi- 
nates Obviously  some  geometrical  information  must  be  put  manually  into  the  computer 
before  it  can  begin  to  discretize  the  cranium  and  its  contents.  However,  from  a minimum  of 
user-specified  input  data,  this  method  resulted  in  an  automatically  computed  and  optimally 
uniform  discretization  of  the  skull-brain  continuum.  There  may  be  some  loss  of  capability 
for  capturing  detailed  geometry,  but  the  overall  ease  and  expediency  of  the  method  far  out- 
weighs that  disadvantage.  I lowever,  the  big  disadvantage  of  the  automatic  scheme  approach 
was  the  additional  effort  required  to  develop  the  necessary  methodology  and  computer  pro- 
gram, which  did  not  exist  at  the  time. 

I he  automatic  skull-brain  mesh  generator  relies  on  Laplace's  equation  in  three  dimen- 
sions, 


V24/(X,Y,Z)  0 (5-1) 

to  uniformly  discretize  the  intracranial  contents.  I his  classical  form  is  applied  to  the  l-J-K 
integer  space  of  the  conceptualized  brain  discretization  shown  in  f igure  5-7.  The  function 
T is  replaced  by  the  unknown  rectangular  Cartesian  coordinates  X,  Y , and  Z,  which  describe 


85 


I 


4 

the  skull-brain  system.  These  functions  are  now  written  in  terms  of  the  independent  integer 
variables  I,  J,  and  K forming  three  independent  Laplace  equations, 


V2X(I,J,K)  = 0 

(5-2) 

V2Y(1,J,K)  = 0 

(5-3) 

V2Z(1,J,K)  = 0 

(5-4) 

These  three  partial  differential  equations  arc  then  cast  into  finite  difference  form,  and 
each  is  solved  separately  using  standard  iterative  techniques  as  described,  for  example,  by 
Smith59. 

As  an  illustration,  the  X -coordinate  function  of  Equation  5-2  is  found  by  solving  a sys- 
tem of  difference  equations  comprised  of  one  equation  for  each  internal  node.  This  results 
in  a system  of  125  equations  (5  cubed)  for  the  example  of  Figure  5-7.  If  a central  differ- 
ence scheme  is  employed  using  the  coordinates  shown  in  Figure  5-8,  a typical  equation  in 
the  system  can  be  written  as 


((N) 

LJ.K 


1 />-l)  V(N)  (N-l)  (N)  (N-l)  (N)  \ 

6 \ I+1,J,K  1-1  J,K  IJ+1.K  IJ-l.K  I,J,K+1  LJ,K-y 


(5-5) 


The  superscripts  (N,N-1,  etc.)  represent  the  iteration  number  during  the  convergence  pro- 
cess. Typically  this  number  can  be  between  30  and  40,  depending  upon  the  prescribed  con- 
vergence tolerance  e;  i.e.. 


(5-6) 


It  is  obvious  from  Equation  5-5  that  the  Laplace  method  results  in  a simple  averaging  for- 
mula for  computing  the  coordinates  of  internal  nodes.  This  results  in  an  optimally  uniform 
brain  discretization.  It  is  worth  noting  that  Gallagher14  has  observed  that  mesh  uniformity 
is  desirable  for  accuracy  in  the  finite  element  solution  process.  This  is  particularly  true  with 
some  incompatible  elements. 


59(,  I)  Smith. 


Numerical  solution  of  partial  differential  equations. 


London,  Oxford  University  Press,  1965. 
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During  the  iteration  proeess  the  neighboring  nodes  surrounding  an  internal  node  near 
the  surface  will  periodically  lay  on  the  cube’s  surface.  The  coordinates  of  such  nodes  are 
prescribed  data  tor  the  Laplace  scheme  and  are  fixed  and  not  subject  to  change  during  iter- 
ation. Put  another  way,  these  data  form  the  boundary  conditions  for  the  solution  of  liqua- 
tion 5-1 . These  data  are  taken  from  actual  skulls  or  plastic  skull  replicas  as  shown  in  Figure 
5-9.  Because  accurate  geometry  of  both  the  skull  and  brain  surface  is  sought,  the  prescribed 
geometrical  data  are  recorded  for  the  internal  skull  surface  as  opposed  to  the  external  sur- 
face. Along  with  the  X-,  Y-,  and  Z-coordinate  values  of  a typical  surface  node,  the  local 
skull  bone  thickness  is  also  recorded.  Graphical  results  of  the  iteration  are  shown  in  Figure 
5-10. 

I he  skull-brain  mesh  generator,  besides  discretizing  the  brain  with  the  Laplace  scheme, 
discretizes  the  skull  bone  into  three  layers  representing  the  inner  and  outer  table  bone  and 
the  middle  layer  of  diploe.  I'his  is  accomplished  directly  (without  iteration)  by  computing 
the  outward  normal  at  each  prescribed  node,  assigning  its  length  to  be  the  skull  thickness 
and  then  appropriately  dividing  the  length  into  three  segments.  Everywhere  the  inner  and 
outer  table  thickness  is  one-fourth  the  skull  thickness,  and  the  diploe  is  half  the  skull  thick- 
ness. A fourth  layer  is  also  provided  at  the  same  time  (using  the  inward  normal)  and  simu- 
lates the  subarachnoid  space  between  the  brain  and  skull.  Its  thickness  is  arbitrarily  selected 
as  being  everywhere  equal  to  the  table  bone  thickness.  The  conceptualized  skull-brain  dis- 
cretization is  shown  in  Figure  5-11. 

Typical  skull-brain  discretizations  that  result  from  the  skull-brain  mesh  generator  are 
depicted  in  f igure  5-12.  The  coarse  facial  bone  discretization  was  constructed  manually 
subsequent  to  the  automatic  generation  of  the  mesh.  Scale  factors  arc  provided  to  alter  the 
shape  or  size  of  an  existing  discretization  (within  certain  limitations),  thus  avoiding  the  nec- 
essity' for  having  to  re-gencratc  the  mesh  when  studying  the  influence  of  skull  size  and 
shape.  If  warranted,  more  detail  in  the  brain  can  be  provided  by  prescribing  the  coordinates 
of  internal  nodes  to  coincide  with  intracranial  structures  or  more  detailed  brain  surface  shape. 


88 


(b)  Three-quarter  view. 

ligure  5-9.  Source  of  Mesh  Generator  Input  Data. 
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figure  5-1  1 Conceptualized  Skull-Brain  Discretization. 
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iidsagittal  plane 


(il)  HIM  external  skull  shape  profile  continuous  tone  computer 
graphics  display. 
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(e)  HIM  external  skull  and  brain,  right  rear  view  continuous  tone 
computer  graphics  display. 


I igtir<  ' i ree  Dimensional  Skull  Brain  Discretization 
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Kornhauscr6"  has  suggested  that  with  proper  validation  the  advanced  computer  or 
mathematical  model  approach  is  the  ultimate  quantitative  prediction  method  in  biodynamic 
modeling.  1'here  is,  however,  very  little  information  on  past  what  “proper  validation”  en- 
tails. Verification  of  biomechanical  models  is  more  complex  than  verification  of  inanimate 
system  models.  I his  is  because  it  must  inevitably  deal  with  experimentally  obtained  in  vivo 
or  in  vitro  data,  where  the  experimental  technique  and  repeatability  of  data  are  adversely 
affected  by  physiological  processes. 

Exact  solutions  are  available  only  for  simple  problems,  and  oftentimes  experimental 
results  are  questionable.  Verification  or  making  sure  the  program  does  what  it  is  intended 
to  do  is  the  responsibility  of  the  originator.  Verification,  however  is  a haphazard  process; 
there  is  no  such  thing  as  absolute  verification.  It  has  been  proposed  by  Griffin61  that  verifi- 
cation be  replaced  by  “level  of  confidence.”  The  level  of  confid  ,,ce  of  a program  would 
encompass  a number  of  things,  including  qualifications  of  the  programmer,  the  complexity 
of  the  model,  determination  of  whether  solution  procedures  are  standard  or  advanced,  de- 
gree of  approximation,  comparison  with  theory  and  experiment,  and  the  length  of  time  the 
program  has  been  in  use.  Ultimately,  however,  it  is  the  user  who  is  the  final  debugger  and 
he  should  not  accept  any  program  as  being  absolutely  verified  but  should  scrutinize  all  the 
data  and  then  feed  back  to  the  program  developer  any  errors  noted. 

I he  linear  head  injury  model  development  was  viewed  as  primarily  art  application  of 
existing  technology  and  consequently  use  of  available  computer  programs  was  indicated  as 
opposed  to  the  development  of  a distinctly  new  code.  The  computer  program  upon  which 
the  HIM  code  is  based  was  developed  by  Taylor26  and  is  termed  the  IT.AP  (f  inite  Element 
Analysis  Program)  code.  This  section  was  primarily  based  on  the  IT.AP  code’s  potential  for 
allowing  the  specification  of  a genera!  purpose  code  of  moderate  size  to  the  unique  problem 
of  head  injury  modeling.  It  was  not  large  and  unwieldy  and  possessed  the  generality  be- 
lieved sufficient  lor  the  task.  Modifications  of  both  an  unforeseen  and  an  anticipated  nature 
could  be  effected  by  delving  into  the  code.  The  element  library  contained  standard  two 
and  three-dimensional  isoparametric  elements,  but  the  IT.AP  code’s  unusual  capability  for 
combining  these  elements  with  linear  viscoelastic  material  behavior  was  decisively  influen- 
tial Many  different  structural  analyses  had  been  successfully  performed  in  the  past  with  the 
I T AP  code,  and  it  was  sufficiently  well-documented.  Since  its  modification  to  the  HIM 
code  form,  there  have  been  close  to  90  static  and  dynamic  head  injury  simulations  per- 
formed Many  of  these  were  developmental  and  of  a debugging  nature.  To  date  there  arc- 
no  known  errors  either  in  the  logic  or  in  the  coding,  f or  the  past  few  years,  the  finite  ele- 
ment method  has  been  established  as  a viable  three-dimensional  analysis  technique.  Two  of 
the  HIM  code’s  features  recognizable  geometry  and  linear  viscoelasticity  arc  well  within 


fit  I 

M Kornhauscr  Biodynamu  modeling  and  scaling  anthropomorphic  dummies,  animals,  and  man,"  paper  presented  at 
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the  state-of-the-art.  Ihe  equations  of  motion  are  linear  and  are  solved  by  standard  solution 
techniques.  The  finite  element  formulations  (element  types)  are  also  standard  eight-node, 
compatible,  isoparametric  elements.  Indeed  these  standard  features  were  initially  selected 
to  enhance  the  level  of  confidence  of  the  HIM  code. 

Correlation  with  experimentally  measured  data  in  primates  has  been  from  the  beginning 
an  objective  in  the  HIM  code  development.  I'he  results  thus  far  have  not  been  completely 
successful.  The  primary  difficulty  is  that  the  experimentally  measured  parameters  are  basi- 
callv  different  in  kind  that  the  computed  parameters.  Nevertheless,  all  of  the  efforts  to 
verify  the  HIM  code  by  correlation  with  experimental  data  arc  presented  and  discussed  in 
this  chapter.  All  of  the  experimental  data  were  supplied  by  the  Highway  Safety  Research 
Institute  (HSRI)  at  the  University  of  Michigan  under  contract  DOT-HS-03 1 -3-749  with  the 
Department  of  Transportation,  National  Highway  Traffic  Safety  Administration. 

Static  Skull  Hone  Structural  Validation 

Static  validation  consists  of  comparisons  between  two  kinds  of  measured  and  computed 
data;  load-deflection  data  and  load-strain  data  for  an  empty  rhesus  skull. 

A midsagittal  plane  view  of  the  actual  rhesus  skull  is  shown  in  f igure  6-1  The  three 
directions  and  regions  of  loading,  A-P,  I.-R,  and  S-I,  arc  shown  relative  to  the  strain  gage 
rosette  locations  employed  by  HSRI.  In  f igure  6-2,  the  midsagittal  plane  view  of  the  cor- 
responding three-dimensional  finite  element  model  is  shown.  Although  not  shown,  the 
facial  bone  discretization  was  actually  included  in  the  computer  simulation  HSRI  provided 
Chd.  with  the  actual  rhesus  skull  from  which  the  experimental  data  was  obtained  I his  skull 
was  discretized  as  described  previously  for  human  skulls  by  the  automatic  mesh  generator. 

The  numbered  finite  elements  shown  were  those  whose  computed  strain  response  was 
compared  with  the  measured  strain  response  from  the  strain  gage  rosettes.  1 wo  elements, 
468  and  540,  arc  designated  in  the  occipital  bone  region  because  the  corresponding  strain 
gage  rosette  is  on  their  common  border. 

To  account  for  the  variance  among  published  mechanical  property  data  for  skull  bone, 
two  limiting  values  were  used  for  Young’s  modulus.  Thus,  computed  data  are  designated 
either  as  corresponding  to  a “high  range”  modulus  or  a “low  range”  modulus.  I he  high 
value  is  1 .78  x 1 (ft  psi  and  is  taken  from  McElhancy  ct  al6’  directly  for  compact  bone.  The 
low  value  is  0.82  x 10^  psi  and  is  derived  as  follows.  Instead  of  compact  bone  properties, 
it  is  based  on  experimental  composite  skull  bone  properties  which  include  the  porous  diploe 
layer  along  with  the  compact  bone.  Additionally  its  derivation  excludes,  as  being  too  stiff, 
all  mechanical  property  data  obtained  from  embalmed  material,  f urthermore,  its  premise  is 
that  the  simulation  of  three  homogeneous  elements  through  the  thickness  must  behave  on 
the  whole  as  does  the  composite  bone.  Mechanical  properties  for  the  middle  skull  layer 
(diploe*)  are  taken  as  one-tenth  the  values  for  the  inner  and  outer  layers. 


|,nnes  11  Mel  111  alley,  et  al.  "Mechanical  properties  ot  cranial  bone,”  Journal  of  Biomechanics,  vol,  3,  1970,  pp 
495-51 1 
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A comparison  of  computed  and  measured  load-deflection  results  for  the  A-P,  l.-R,  and 
s I directions  is  presented  in  Figure  6-3.  I he  loading  portion  of  the  experimental  curves 
demonstrate  the  linearity  of  actual  skull  deflections  within  the  range  of  applied  load.  ( I he 
loading  rates  varied  from  3 to  9 Ib/scc.)  However,  the  unloading  portion  consistently  ex- 
hibited a nonlinear,  hvsteritic  behavior.  I he  HIM  code  is  linear  and  does  not  account  for 
this  phenomenon.  I he  essence  of  this  data  involves  a comparison  between  measured  and 
computed  slopes  (average  of  the  slopes  for  experimental  data).  These  values  are  a measure 
of  the  skull's  structural  stiffness  in  each  of  three  roughly  orthogonal  directions. 

A summary  of  the  data  is  presented  in  Table  6-1  and  shows  clearly  that  the  finite  cle- 
ment model  is  from  1 .2  to  2.2  times  more  stiff  for  the  lower  modulus  and  from  2.7  to  4.7 
times  more  stiff  for  the  higher  modulus  than  the  actual  skull.  At  first  it  was  believed  that 
this  discrepancy  was  due  to  errors  in  simulated  skull  thickness  values.  Hut  after  carefully 
reconstructing  the  geometrical  simulation,  only  a small  improvement  was  achieved.  It  is 
now  believed  that  the  difficulty  is  w ith  the  particular  finite  element  formulation  used  for 
the  skull  element. 


Table  6-1 . Comparison  of  Experimental  and  Computed  Structural  Stiffness 

for  Rhesus  Skull  No.  1 


Direction 

I °[ 

1 oading 

HSKI  Experimental 
Stiffness 
(Ib/in) 

Computed  Stiffness 
Using  l.ow  E 
(lb/in) 

Computed  Stiffness 
Using  High  E 
(lb/in) 

A-P 

12,600 

15,527 

33,69  3 

l.-R 

3,700 

8,01 1 

17,367 

SI 

3,087 

5,777 

12,470 

I he  element  employed  is  the  eight-node  brick,  and  it  cannot  effectively  capture  large 
trending  stress  gradients.*  In  the  region  of  the  skull  bone  directly  beneath  the  load  the  de- 
lormation  is  primarily  due  ter  bending.  However,  in  regions  only  a small  distance  away  the 
bending  stress  gradients  become  small  and  the  membrane  stress  values  predominate.  Here 
the  element  formulation  is  expected  to  work  well.  But  in  areas  of  large  bending  stress 
gradients,  the  solution  is  to  use  a different  finite  element  formulation;  either  a 20-node 
brick,  or  a nonconforming  formulation.  To  correct  the  problem  an  eight-node,  incompat- 
ible finite  element  has  been  developed  and  checked  by  'Taylor,  et  al.6j  and  has  been  added  to 


'Several  concomitant  cantilever  beam  model  studies  verified  this 

Robert  I laylor.  Peter  I liereslord,  and  l.dward  I Wilson  \ nonconforming  element  for  stress  analysis,"  Inter 
national  Journal  lor  Numerical  Methods  in  I ngmeenng  (submitted  lor  publication,  l‘>75). 
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the  I I Al  codi  element  library  I he  element  has  demonstrated  improved  bending  response 
and  has  been  coded  to  list  strains  on  the  element  surfaces.  Unfortunately  funds  expired 
prior  to  incorporation  of  this  element  into  the  HIM  code. 

Comparisons  of  computed  ami  measured  strains  are  presented  in  figure  6-4  The  data 
are  shown  in  "load  versus  principal  strain"  format.  1 oads  were  experimentally  applied 
quasi  statically  in  15  tests  which  are  numbered  consecutively  from  3 to  17.  Basically,  the 
experimental  data  is  acceptable  and  once  again  demonstrates  the  linearity  of  skull  bone  for 
the  range  ol  loading  employed.  In  some  load  cases  the  spread  in  data  is  larger  than  in 
others.  A trained  eye  may  be  able  to  discount  some  of  it,  but  all  the  data  received  from 
HSR  I is  plotted  in  this  validation  effort.  I he  computed  data  is  shown  for  both  the  high  and 
low  elastic  modulus  h with  a solid  line  and  a dashed  line,  respectively.  I ach  line  is  num- 
j bered  corresponding  to  the  element  from  which  the  computed  principal  strain  w'as  ex- 

tracted. Also  in  each  figure  a Mohr’s  circle  for  strain  is  qualitatively  sketched  for  both  the 
measured  and  computed  strain  states,  f rom  figure  to  figure  the  relative  circle  si/e  describes 
qualitatively  the  relative  strain  magnitudes  occurring  from  test  to  test.  I he  solid  semicircle 
represents  the  measured  values. 

I he  primary  consideration  in  looking  over  the  strain  comparisons  is  that  the  measured 
values  are  surface  strains  and  the  computed  values  are  subsurface  strains  taken  from  an  ele- 
ment in  the  vicinity  of  the  strain  gage.  Within  the  element  the  computed  strain  is  every- 
where constant,  but  it  portends  to  be  most  accurate  at  the  element’s  centroid.  I he  result 
j is  that  comparisons  are  being  attempted  of  measured  and  computed  strains  at  points  which 

differ  in  location ; i.e.,  the  gage  point  and  the  element  centroid.  Depending  on  the  strain 
variation  in  the  neighborhood  of  these  two  points,  this  difference  may  or  may  not  be  sig- 
nificant Generally,  however,  it  is  known  that  the  shear  strain  at  subsurface  position  is 
nonzero  (tends  towards  a maximum  at  the  midthickness  point)  whereas  on  the  free  surface, 
the  two  shear  strain  components  must  always  be  zero  in  the  absence  of  a traction  force. 
Since  shear  strain  is  equivalent  to  the  Mohr  circle  radius,  the  computed  Mohr  circles  (for 
subsurface  strains)  are  always  expected  to  be  larger  than  the  measured  Mohr  circles  and  the 
results  do  bear  this  out. 

In  spite  of  the  above  limitations,  at  least  qualitative  evaluation  can  be  made  of  the  com- 
parison presented  in  the  nine  parts  ol  figure  6-4.  Ihe  results  range  from  acceptable  corre- 
lation to  no  correlation  at  all.  I bis  qualitative  evaluation  is  summarized  in  Table  6-2. 

figures  6 4(d)  and  (e)  show  the  computed  strains  bracket  the  experimental  strains 
nicely,  f his  is  true  to  a somewhat  lesser  extent  in  f igures  64(c)  and  (i)  if  we  choose  to  key 
on  element  468  instead  of  540.  In  all  cases,  element  468  seems  to  correlate  better  than 
element  540.  I hus,  in  these  four  cases  if  an  average  Young’s  modulus  were  employed  the 
( computed  strains  would  correlate  acceptably  well  with  average  measured  strains.  This  could 

also  be  true  of  the  results  in  I igure  64(a)  although  the  experimental  data  spread  simply 
preclude-  any  meaningful  evaluation.  I lie  results  show  correlation  ranging  from  poor  to 
none  at  all  in  figures  64(b),  (f),  (g),  and  <h).  Results  of  figures  64(b)  and  (h)  demonstrate 
a qualitative  agreement.  Again  assuming  an  average  modulus,  their  error  ranges  from  50  to 
1°0%  However,  the  results  shown  in  f igures  64(f)  and  (g)  show  the  model  to  be  predicting 
much  more  strain  (more  than  100%)  than  what  is  being  measured. 
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Dynamic  Kxperimcntal  Validation 

I here  are  two  kinds  of  data  and  two  different  configurations  involved  in  the  dynamic 
validation  phase  which  immediately  followed  the  static  validation  effort.  Dynamic  cranial 
bone  strains  in  both  a drv  rhesus  skull  and  a live  rhesus  were  to  be  used  as  correlative  data. 
These  strains  were  induced  by  impacts  in  a controlled  laboratory  environment.  More  im- 
portantly, dynamic  intracranial  pressure  data  recorded  for  the  live  animal  impacts  were  to 
be  used  to  correlate  measured  and  predicted  brain  response. 

Dynamic  Skull  Hone  Response.  In  the  first  configuration  the  empty  rhesus  skull  was 
anchored  to  a laboratory  test  frame  with  a wall-anchor  bolt  assembly  inserted  through  the 
skull's  foramen  magnum  so  that  the  skull  base  was  immovable  against  various  dynamic  im- 
pacts. lour  strain  gage  rosettes  were  located  on  the  frontal,  the  two  parietal,  and  the  occi- 
pital cranial  bones,  t he  HSRl  principal  strain  data  shown  in  figure  6-5  were  recorded  for 
a 3 .8  ms  impact  on  the  superior  surface  in  the  S-l  direction.  The  peak  load  was  63  pounds 
applied  over  an  area  ol  1/2  sq.  in.  The  strain  magnitudes  are  typical  of  this  series  of  exper- 
imental tests.  1 he  first  impression  of  these  measured  strains  were  that  they  appeared  too 
low.  They  are  on  the  whole  about  five  times  smaller  than  the  strains  recorded  during  the 
static  validation  tests  at  equivalent  peak  load  levels.  This  difference  can  be  explained  from 
an  examination  of  the  geometry  of  the  applied  load  and  boundary  condition  at  the  skull 
base. 

In  the  static  load  cases  the  applied  S-l  loading  was  equilibrating,  causing  the  rhesus  skull 
to  be  “pinched'1  in  the  S I direction  (see  figure  6-1 ).  In  the  dynamic  case  the  applied  S-l 
load  vector  does  not  pass  through  the  base  support  but  rather  forward  of  it,  and  instead  of 
pinching  the  skull  it  bends  the  skull  as  in  a cantilever.  The  configuration  of  the  correspond- 
ing HIM  code  simulation  illustrates  this  in  f igure  6-6.  I his  causes  rigid  body  motion  of  the 
skull  about  its  support,  but  motion  not  present  in  the  static  tests.  Thus,  the  dynamic  strain 
levels  recorded  may  well  be  below  those  recorded  for  the  static  case  as  mentioned  above. 

If  the  test  setup  had  been  arranged  so  that  the  dynamic  load  vector  passed  through  the 
support,  subsequent  simulation  of  the  test  would  have  gone  more  smoothly.  Review  of 
I1SRI  high  speed  film  revelacd  that  the  wall-anchor  bolt  assembly  deflected  elastically  under 
the  impact  of  the  controlled  blow.  This  meant  that  the  corresponding  HIM  code  simulation 
had  to  account  for  this  motion  instead  of  merely  constraining  the  nodal  points  at  the  fora- 
men magnum  against  movement.  Several  load-deflection  tests  of  the  wall-anchor  bolt  as- 
sembly had  to  be  conducted  at  ( Id.  to  determine  its  structural  properties.  These  proper- 
ties were  then  translated  into  an  equivalent  wall-anchor  bolt  finite  element  which  was  added 
to  the  base  of  the  simulted  skull.  This  process  is,  at  best,  only  roughly  approximate.  Never- 
theless, the  computed  strains  for  the  simulated  S I impact  are  shown  in  f igure  6-7.  I he 
upper  curve  represents  the  computed  strain  history  for  complete  restraint  of  the  skull  base, 
and  the  lower  curve  represents  the  strain  history  for  an  elastically  restrained  skull  base.  Ob- 
viously the  restraint  greatly  influences  the  predicted  strain.  Although  many  more  computer 
simulations  were  conducted  in  this  validative  effort,  the  two  curves  represent  the  limits  of 
strain  response  predicted. 
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l-igurc  6-5  I-  xperimcntai  Dynamic  Strain  in  a Klu-sus  Skull 
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l igure  6-7  Computed  Dynamic  Strain  in  Kiicsus  Skull  Parietal  Hon 


I lie  primary  observation  when  comparing  the  measured  anil  computed  strains  is  that  the 
computed  strain  is  at  least  three  times  greater  in  magnitude,  ((.'are  was  taken  that  the  com- 
puted strain  was  sampled  at  a location  near  the  corresponding  strain  gage  rosette  on  the 
rhesus  parietal  bone.)  Comparisons  of  computed  and  measured  dynamic  strains  elsewhere 
on  the  rhesus  skull  consistently  revealed  that  the  model  predicts  strains  by  a factor  between 
three  and  one  order  of  magnitude  greater  than  measured  strains.  The  discrepancy  was 
ascribed  to  the  eight-node  isoparametric  element  whose  stiffness,  as  discussed  in  the  static 
validation,  did  not  adequately  capture  the  skull  bending  response.  Because  of  the  inherent 
difficulties  in  attempting  to  match  finite  element  model  strains  with  those  of  the  actual 
skull,  it  was  expected  that  correlation  would  be  difficult.  I he  finite  element  model,  how- 
ever, should  be  capable,  in  this  instance,  of  better  correlation  and  on  that  basis  requires 
improvement. 

Dynamic  Intracranial  Pressures.  A series  of  experimental  validation  tests  was  conducted 
by  HSR1  with  live  rhesus  monkeys  in  an  anesthetized  condition.  In  the  case  of  A-P  (rear) 
impacts,  two  specially  prepared  pressure  transducers  were  surgically  placed  in  the  skull  to 
monitor  intracranial  pressure.  One  was  located  in  the  left  parietal  bone  near  the  point  of 
impact,  and  the  other  was  in  the  forward  right  corner  of  the  frontal  bone  over  the  tip  of  the 
right  frontal  lobe  of  the  brain.  A schematic  of  the  pressure  transducer  system  is  presented 
in  f igure  6-8.  I lie  transducer  sensor  is  not  located  in  a subdural  position;  instead,  it  is  care- 
fully placed  epidurally  but  very  close  to  the  dura.  Approximately  one  eve  drop  of  silica  gel 
(ills  the  notch  in  the  skull  and  is  intended  to  provide  a pressure  transmitting  medium.  In  this 
way  the  pressure  measured  is  said  to  be  equivalent  to  the  intracranial  pressure  even  though 
the  sensor  is  not  in  direct  contact  w ith  any  intracranial  material. 

A typical  measured  pressure  history  is  shown  in  f igure  6-9.  In  this  case  the  data  are 
from  the  front  pressure  transducer,  and  the  applied  load  was  from  the  rear  so  that  the  data 
are  said  to  be  contrecoup  pressure,  flic  applied  peak  load  magnitude  in  HSR  I run  number 
o<)3  was  800  pounds  and  the  duration  ol  loading  was  2 ms  which  amounts  to  fairly  rigid 
impact  I he  initial  peak  pressure  is  about  1 3 psi  tension  and  occurs  at  about  1 ms  and  the 
next  peak  is  19  psi  compression  occurring  just  after  2 ms.  The  second  peak  is  believed  due 
to  the  influence  of  the  rhesus  neck  structure  in  reversing  the  forward  acceleration  due  ini- 
tially to  the  applied  load.  Since  the  primate  is  heavily  anesthetized,  this  may  be  neither  a 
voluntary  nor  an  involuntary  neck  reaction  but  just  from  the  physical  presence  of  the  neck 
and  thoracic  inertia.  I lie  first  peak  is  the  classical  tensile  pressure  at  the  contrecoup  loca- 
tion More  importance  is  placed  on  it  because  it  is  popularly  thought  that  tensile  pressures 
arc  more  conducive  to  brain  injury. 

Several  associated  computer  runs  with  the  HIM  code  were  made  in  an  attempt  to  corre- 
late the  computed  response  with  the  above  measured  intracranial  response.  In  each  succes- 
sive computer  run  adjustments  to  the  model  were  made  which  were  believed  to  be  improve- 
ments in  either  the  model  or  the  simulation  and  which,  it  was  hoped,  would  bring  the  com- 
puted pressure  into  closer  agreement  with  the  measured  pressure. 

I he  midsagittal  plane  discretization  ot  the  simulated  rhesus  skull  is  shown  in  f igure 
6 10  Neither  the  scalp  and  other  soli  tissue  nor  the  mandible  are  included  directly  in  the  dis 
erctization  I lie  combined  mass  of  these  components  was,  however,  subsequently  com- 
puted and  added  to  the  model  by  increasing  the  density  of  the  outer  layer  of  skull  bone 
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Note  1 Not  to  scale. 

2.  Notch  is  Irom  7 to  10  mm  in  diameter 

3.  Sensing  element  tip  is  about  2 or  3 mm  off  the  dura  surface 
4 Only  an  "eve  drop  of  silica  gel  is  inserted  into  notch. 
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figure  6 K Schematic  of  I IS  HI  Pressure  Transducer  and  Mount. 


material  by  an  appropriate  uniform  amount.  The  facial  bone  structure  of  the  rhesus  skull 
is  very  complex  and  extensive  when  compared  to  the  human  skull.  Although  much  effort 
was  expended  in  an  attempt  to  include  the  mandible  in  the  rhesus  discretization  it  proved 
too  difficult  (unlike  the  human  model  in  which  the  entire  facial  bone  complex  was 
included).  The  other  facial  bones  were  difficult  and  their  discretization  (shown  in  figure 
6-10)  was  not  achieved  easily.  Perhaps  a more  important  consideration  in  evaluating  the  live 
rhesus  monkey  simulation  configuration  is  the  unknown  boundary  condition  at  the  base  of 
the  skull.  As  stated  above,  the  neck  is  believed  to  influence  the  rhesus  head  acceleration 
(actually  reversing  it)  subsequent  to  the  impact.  Since  the  corresponding  model  does  not 
possess  a neck  structure,  the  neck  reaction  forcing  function  on  the  skull  would  normally  be 
considered  as  necessary  input  information  as  is  the  simulated  impact  force.  If  this  informa- 
tion were  available  it  could  easily  be  input  to  the  model  nodes  at  the  base  of  the  skull  but 
it  is  not  known  nor  is  it  known  how  it  can  be  measured.  Therefore,  the  only  recourse  is  to 
assume  some  reasonable  displacement  boundary  condition  for  the  nodes.  The  choices  for 
this  condition  in  the  present  HIM  code  version  are  completely  unrestrained  nodes,  com- 
pletely restrained  nodes,  or  constrained  nodes  that  slide  along  a prescribed  flat  plane. 

(These  cases  are  evaluated  and  discussed  in  conjunction  with  the  human  model  in  the  next 
chapter.)  The  best  course  ot  action  appears  to  be  to  assume  a completely  unrestrained 
(free)  condition  for  these  nodes.  At  least  the  initial  peak  (tensile  peak)  in  the  measured 
intracranial  pressure-history  should  be  unaffected  by  the  neck  reaction.  Indeed  the  IISRI 
experimental  configuration  was  designed  to  create  a ‘‘free  body”  condition  for  the  live 
rhesus  head  at  impact  by  aligning  the  head  in  front  of  the  impactor  with  light  supportive 
surgical  thread  attached  to  the  primate's  ears.  Thus,  it  is  believe  to  be  better  in  the  corres- 
ponding simulation  to  model  the  rhesus  skull  with  no  restraint  at  the  skull  base.  Computed 
intracranial  pressure  data  was  extracted  and  examined  at  element  109,  which  corresponds 
closely  to  the  HSK1  transducer  location,  and  also  along  an  anterior-posterior  (A-P)  plane  as 
shown.  I he  anterior  point  (A),  intermediate  point  (I),  ami  posterior  point  (P)  are  in  the 
midsagittal  plane.  I lenient  109  is  not  in  the  midsagittal  plane. 

I he  simulated  impact  force-history  is  shown  in  figure  6-1 1 and  was  uniformly  distri- 
buted over  an  area  of  approximately  1/2  sq  in  It  simulates  the  short  duration  blow  deliv- 
ered symmetrically  about  the  midsagittal  plane  with  a hard  impactor. 

Results  of  the  HIM  code  simulations  of  the  live  rhesus  monkey  test  are  shown  in  f igures 
612(a)  through  <f)  An  integration  time  step  size  of  1 30  /as  was  employed  in  all  the  compu- 
tations except  as  noted  I Ins  series  of  six  computer  runs  consitute  the  more  important  at- 
tempts to  compute  intracranial  pressures  which,  it  was  hoped,  would  correlate  with  those 
pressures  measured  (f  igure  6-9).  More  specifically  the  objective  of  the  computation  was  to 
agree  as  closely  as  possible  with  the  peak  tension  pressure  of  1 3 psi  observed  by  IISRI. 

In  the  first  simulation,  f igure  6 12(a),  the  rhesus  skull  model  contained  no  opening 
through  the  cranial  base  where  the  foramen  magnum  is  normally  found  and  through  which 
the  brainstem  passes.  I his  was  the  existing  configuration  of  the  model  at  that  time. 

Though  element  109  produced  a relatively  high  peak  tensile  pressure  of  3 5 -psi  tension,  these 
results  were  initially  encouraging  because  it  was  believed  that  the  intended  subsequent  in- 
clusion of  the  foramen  magnum  would  uniformly  lower  the  computed  pressure.  Also  it 
was  very  uncommon  to  find  any  head  injury  model  predicting  peak  intracranial  pressures 
below  100  psi  at  that  time.16 
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(I))  foramen  magnum  included. 


(c)  Compressible  subarachnoid  space. 
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!n  the  remaining  five  simulations  the  specification  of  bone  material  properties  for  the 
finite  elements  in  the  region  of  the  foramen  magnum  (see  f igure  6-10)  was  removed  and  re 
plaeed  as  follows.  I he  outer  layer  of  elements  were  assigned  to  be  sold  and  elastic  in  the 
manner  of  a membrane  sealing  over  the  foramen  magnum  aperture.  I he  nest  three  inner 
most  layers  were  assigned  the  mechanical  properties  of  the  brain  material  and  were  thought 
of  as  simulating  a rather  truncated  brainstem.  In  the  absence  of  a neck  model,  it  was  hoped 
that  this  would  provide  a more  realistic  approximation  to  the  region  of  the  foramen  mag 
mini  and  brainstem  while  still  providing  containment  (necessary  for  numerical  stability)  for 
the  simulated  intracranial  contents. 

I he  results  in  figure  6-1  2(b)  were  very  surprising,  however,  because  while  the  computed 
pressure  at  point  lJ  decreased,  the  tensile  pressure  extracted  from  element  109  increased 
front  35  psi  to  67  psi  A general  decrease  in  overall  intracranial  pressure  did  not  occur  as 
was  expected,  but  instead  a shif  t in  the  distribution  of  the  pressure  materialized.  1'his  was 
unfortunate  from  the  standpoint  of  achieving  correlation  of  the  1 3 psi  peak  tensile  pressure 
in  the  { egion  of  element  109  But  tor  reasons  which  will  be  explained  more  fully  in  a sub 
sequent  discussion  on  statistics  of  observed  head  injury  distribution  a return  to  the  “no 
foramen  magnum"  configuration  could  not  be  justified  Also  in  this  simulation  the  applied 
load  was  improved  by  including  both  I | and  f ^ exactly  as  shown  in  f igure  6-1 1 whereas 
the  previous  simulated  load  included  onl\  I j.  Comparison  w ith  the  measured  impact  force 
history  showed  that  the  combination  of  I j and  I ? gave  a better  approximation  of  the  ac- 
tual load  I his  change  would  explain  the  difference  between  f igure  6 1 2(a)  and  (b)  in  the 
“trailing  off”  portion  of  the  computed  intracranial  pressures.  Something  still  had  to  be 
done  w ith  the  model  to  reduce  the  predicted  pressures,  and  there  remained  a modification 
that  had  been  intended  but  yet  w;as  untried. 

Since  the  brain  is  not  continuous  with  the  internal  cranial  wall  but  is  supported  and 
tethered  at  the  skull/brain  interface,  a method  was  devised  to  approximate  this  suspension. 

I he  single  layer  of  elements  which  simulate  the  subarachnoid  space  was  provided  in  the 
basic  I II.M  code  configuration  for  this  eventuality.  Up  to  this  point,  the  mechanical  prop- 
erty specification  for  these  elements  has  been  identical  to  that  of  the  brain  elements,  f his, 
combined  with  the  displacement  continuity  imposed  by  the  finite  element  method,  effects 
a continuous  support.  1 lowcvcr,  by  specifying  the  material  for  the  subarachnoid  elements 
(see  f igure  6-10)  as  compressible,  a discontinuous  support  between  the  brain  anil  skull  can 
be  approximated.  As  it  happens  the  net  result  should  also  be  to  reduce  the  deformation 
gradients  and,  hence,  the  intracranial  pressures  that  develop  in  the  brain. 

I he  results  of  characterizing  the  subarachnoid  space  as  compressible  (bulk  modulus  of 
3,050  psi)  are  shown  in  f igure  6 12(c).  I his  time,  as  expected,  the  overall  computed  intra- 
cranial pressure  was  reduced  Pressures  at  all  points  were  lower  than  they  were  prior  to  the 
compressibility  specification  in  the  subarachnoid  space  elements.  I he  peak  tensile  pressure 
in  element  109  dropped  front  67  to  55  psi.  I his  was  still  high  compared  to  the  experi- 
mental value  of  1 3 psi  but  nevertheless  was  considered  an  improvement  Strains  in  the 
brain,  though  not  shown,  were  reduced  correspondingly. 

Varying  the  compressible  properties  of  the  subarachnoid  space  is  an  attempt  to  deal 
with  an  intracranial  area  for  which  there  exists  no  material  propertv  data  and  in  which  the 
stress  strain  behavior  is  of  no  interest  so  long  as  discontinuity  of  displacements  is  simulated. 


I here  may  never  be  any  material  property  data  because  of  the  discontinuous  nature  of  the 
subarachnoid  space.  One  recourse  is  to  effect  engineering  approximations  with  regard  to  the 
mechanical  behavior  of  this  area. 

In  the  course  of  continued  cheeking  of  the  simulation  for  accuracy  it  was  thought  that 
the  mass  of  the  connecting  head  tissue  and  the  mandible,  which  was  neglected  in  the  model, 
might  be  more  significant  than  at  first  realized.  An  estimate  of  this  additional  mass  was 
made  for  the  specific  primate  employed  in  IISK1  test  003.  Distribution  of  this  additional 
mass,  when  added  to  the  model,  could  not  be  accomplished  accurately  since,  as  described 
earlier,  no  discretization  of  the  rhesus  mandible  was  included.  But  the  mass  density  of  the 
outer  layer  of  elements  tor  the  cranial  bone  was  increased  by  a uniform  amount,  making 
the  total  mass  of  the  model  more  closely  equal  to  that  for  the  actual  animal. 

To  demonstrate  the  influence  of  this  additional  mass,  the  rhesus  skull  displacement  data 
are  presented  in  figure  6-13.  1 his  information  was  also  derived  to  aid  in  improving  the 
simulation  of  the  live  rhesus  experimental  tests.  If  a reasonable  correlation  between  mea- 
sured and  computed  intracranial  pressure  is  to  occur,  at  least  the  measured  anil  computed 
displacement  histories  of  the  rhesus  skull  should  agree.  Acceleration  histories  would  per- 
haps have  been  more  convincing  parameters,  but  this  information  was  not  of  sufficient 
quality  to  be  usable  at  the  time  of  the  validation  effort  herein  described.  High  speed  film 
proved  more  useful  in  determining  displacements  than  did  the  double  integration  of 
measured  accelerometer  data.  I he  “particle  mechanics  solution”  assumes  perfect  transla- 
tion which  was  not  the  ease  as  some  rotation  was  present  in  the  horizontal  displacement  of 
the  frontal  bone.  The  displacement  history  curve  associated  with  the  corrected  or  added 
mass  is  more  parallel  with  the  experimental  data  than  any  of  the  other  curves.  This  means 
that  the  computed  velocity  agrees  more  closely  with  the  measured  velocity  as  a result  of  the 
added  mass.  1 his  is  not  necessarily  true  of  the  acceleration,  however,  because  that  quantity 
depends  on  the  curvature  and  not  the  slopes  of  the  curves,  and  curvature  cannot  be  as 
easily  observed  in  this  data.  It  does  appear  as  though  the  added  mass  has  retarded,  by  too 
much,  the  displacement  of  the  frontal  bone. 

A distinct  reduction  in  overall  intracranial  pressure  computed  by  the  HIM  code  as  a re- 
sult of  the  added  mass  can  be  seen  in  f igure  6-1 2(d ).  The  peak  tension  stress  in  element 
109  has  been  lowered  from  55  to  43  psi.  I hese  results  demonstrate  the  importance  of  ade 
quately  simulating  the  mass  of  the  rhesus  skull  when  computing  the  intracranial  pressure. 

At  this  point  there  were  few  obvious  additional  modifications  to  the  simulation  and 
model  which  could  reasonably  be  proposed  without  revamping  the  entire  model  < UvviousK 
the  measured  pressures  were  still  well  below  those  which  were  being  computed  Since  time 
did  not  allow  any  extensive  modification  to  the  model,  two  more  simulations  were  con 
ducted  which  required  only  minor  changes, 

I he  mechanical  properties  of  the  subarachnoid  space  being  quite  unknown  could  be 
guessed  at  with  relative  freedom.  If  a reduction  of  the  skull  modulus  from  305.000  to 
3,050  psi  provided  an  improvement,  then  perhaps  a further  reduction  to  305  psi,  for  ex 
ample,  would  also  be  beneficial.  I'he  computed  intracranial  pressures  presented  in  figure 
6 12(e)  are  for  this  case.  An  improvement  was  indeed  realized  although  not  as  significant 
as  before.  The  computed  peak  tensile  pressure  in  element  109,  corresponding  to  the  pres 
sure  transducer,  was  reduced  from  43  to  37  psi.  I here  is  a limit  imposed  by  numerical 
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stability  to  the  amount  by  which  the  bulk  modulus  can  be  reduced  or  to  the  specified  com- 
pressibility tor  the  layer  of  elements  designated  as  the  subarachnoid  space.  I he  latter  bulk 
modulus  value  of  305  psi  is  straining  that  limit.  Any  further  reduction  would  risk  compu- 
tation of  incompatible  displacements  in  the  subarachnoid  region  with  the  present  model. 

I he  final  attempt  at  correlating  with  the  measured  intracranial  pressures  was  made  by 
testing  the  effect  of  the  prescribed  integration  time  step  size.  To  this  point  1 30  /as  was  the 
time  step  employed  in  all  simulation  and  this  value  might  appear  coarse  when  compared 
with  the  rise  time  associated  with  the  input  load  function  (Figure  6-1  1).  I hereforc.  the  pre- 
scribed time  step  was  halved  to  a value  of  65  /as  which  is  on  the  order  of  one-tenth  of  the 
rise  time  possessed  by  the  impact  load.  I his  is  ordinarily  deemed  adequate,  and  anywhere 
from  one-eighth  to  one-tenth  can  be  used  as  a rule  of  thumb.  1 he  best  method,  however, 
is  to  test  the  effect  on  response  of  various  prescribed  time  steps.  1 he  results,  in  this  case, 
arc  presented  in  figure  612(f),  exhibit  very  little  difference  in  peak  response,  and  show 
that  no  real  benefit  can  be  expected  by  undergoing  the  additional  expense  of  a smaller  inte- 
gration time  step  size. 

A summary  of  thedllM  code  simulations  of  the  live  rhesus  impact  experiment  which 
have  been  discussed  aou.  e is  presented  in  fable  6-3.  Here  it  can  be  seen  that  the  best  effort 
to  correlate  computed  anti  measured  intracranial  pressures  results  in  the  HIM  code  predict- 
ing pi  tk  intracranial  pressures  which  are  about  three  times  higher  than  those  which  were 
measured  in  I1SKI  test  number  003.  It  can  be  seen  that  the  predicted  pressure  gradient  in 
the  A F direction  is  substantial.  It  is  about  65  psi  in  the  short  space  of  approximately  3 
inches.  More  detailed  computed  intracranial  pressure  data  for  the  midsagittal  plane  result- 
ing from  HIM  code  simulation  number  5 (see  Table  6-3)  is  presented  in  f igure  6 14  If  at  a 
later  time  the  position  which  was  selected  for  correlation  with  the  HSR1  transducer  is  found 
to  be  inaccurate,  this  data  can  serve  to  adjust  the  above  validation  results.  As  evidenced  by 
the  pressure  gradient  size,  it  would  not  require  much  of  an  error  in  this  position  to  make  a 
significant  difference  in  these  results. 

Qualitative  Validation 

There  arc  some  characteristic  traits  of  typical  HIM  code  predictions  which  have  sur 
faced  during  the  development  process  and  which  are  consistent  with  general  and  experi- 
mental observation  or  head  injury  statistics.  As  a result  they  are  presented  here  to  enhance 
the  level  of  confidence  in  the  performance  of  the  HIM  code  in  a qualitative  sense. 

It  is  generally  agreed  by  engineers  and  clinicians  alike  that  contrccoup  brain  injury  is  an 
important  category  of  head  injury.  When  the  initial  development  of  the  HIM  code  began 
there  was  no  guarantee  that  this  important  mechanism  would  manifest  itself  in  a finite  ele- 
ment model  5 ct  from  the  preliminary  two-dimensional  models  to  the  present  three- 
dimensional  models,  significant  tensile  intracranial  pressure  has  been  consistently  predicted 
in  that  portion  of  the  brain  opposite  to  the  applied  impact  load.  At  the  same  time,  com- 
pressive intracranial  pressure  in  the  brain  adjacent  to  the  applied  load  (coup  pressure)  is 
consistently  predicted.  I Ins  phenomenon  has  not  only  been  observed  clinically  and  patho- 
logy illy , but  it  is  a reasonable  physical  mechanism  which  should  be  expected  from  a prop- 
erly simulated  encapsulated,  fluidlike  material  I he  model  simulates  this  basic  mechanism 
vtrv  well 
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'able  6-3.  Summary  ot  HIM  (lode  Simulations  of  l ive  Rhesus  I \periments 


Some  interesting  statistics  relating  the  incidence  of  coup  and  contrecoup  head  injuries 
have  been  compiled64  It  is  believed  that  this  information  is  germane  to  the  development 
ot  head  injury  models  and  that  one  criterion  for  a valid  model  is  that  it  account  for  this 
data.  I he  information  is  presented  in  fable  6-4.  It  is  readily  apparent  from  this  data  that 
the  traumatic  results  caused  by  frontal  and  rear  impacts  are  different.*  The  data  suggest 
that  frontal  impacts  are  more  apt  to  induce  coup  brain  injury  than  contrecoup  injury,  while 
rear  impacts  induce  predominantly  contrecoup  injury.  Side  impacts  are  more  apt  to  pro- 
duce contrecoup  injury  independent  of  impact  location  (left  or  right  side).  Symmetry  of 
the  resulting  injury  exists  for  side  impacts  but  does  not  exist  for  front  and  rear  impacts. 


1 able  6 4.  Distribution  of  Coup  and  Contrecoup  brain  Injuries  with 
Location  of  Impact 


1 .ocation 
of 

impact 

Number 

of 

cases 

Percent  location  of  injuries 

Only  or 
mainly  at 
impact  site 

Only  or 

mainly  opposite 
impact  site 

Equally 

distributed 

frontal 

280 

49 

5 

46 

Rear 

36 

0 

97 

3 

Left  side 

33 

12 

67 

21 

Right  side 

31 

13 

68 

19 

In  I igure  6-15,  the  intracranial  pressure  distribution  in  the  A I*  direction  arc  plotted  for 
simulations  I and  2 ( I able  6-3).  I'hese  data  are  normalized  on  the  highest  computed  abso- 
lute pressure  I he  contrecoup  pressures  or  tensile  pressures  naturally  occur  in  the  anterior 
brain  while  the  coup  pressures  or  compressive  pressures  occur  in  the  posterior  brain.  Ilow- 
cur.  note  that  the  A I*  pressure  distribution  curves  for  the  closed  skull  intersect  the  zero 
pressure  line  at  a position  considerably  anterior  to  where  the  open  skull  curves  intersect. 
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* \i  the  very  least  then  .1  model  should  lie  capable  ol  distinguishing  between  frontal  and  rear  impacts. 
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l igurc  6-15  Influence  of  l oramcn-Magnum/Intracranial-l'rcssurc  Distribution. 
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1 he  net  effect  is  that  in  the  closed  skull  the  contrccoup  pressures  are  much  smaller  than  the 
contrecoup  pressures  tor  the  open  skull  ease,  since  at  any  time  the  two  eases  yield  nearly 
eipial  gradients.*  II  the  vulnerability  of  brain  material  to  tensile  pressure  is  accepted,  the 
pressure  distribution  occurring  in  the  open  skull  is  far  more  likely  to  explain  the  data  of 
I able  6-4  1 his  data  indicates  that,  tor  rear  impacts,  in  97%  of  the  eases  brain  injury  occurs 

only  or  mainly  at  a contrecoup  point,  and  in  only  3%  of  the  eases  injury  equally  distributed. 
In  no  eases  were  the  coup  injuries  tounel  to  occur  in  the  absence  of  contrecoup  injuries.  Vet 
the  closed  skull  simulation  would  predict,  in  absolute  magnitude,  coup  pressures  much 
higher  than  contrecoup  pressures.  While  it  may  be  intuitively  obvious  to  some,  this  data 
concludes  that  a closed  skull  configuration  for  head  injury  simulations  must  be  rejected. 

Another  possible  explanation  of  the  data  in  Table  64  would  include  the  differences  in 
the  bony  irregularities  ot  the  anterior  and  posterior  internal  skull  surface.  Since  the  anter- 
ior skull  contains  more  sharp  protrusions,  the  interface  conditions  between  the  cranial  wall 
and  brain  and  attending  hazards  will  be  more  severe  under  dynamic  conditions.  This  is  un- 
doubtedly true,  but  the  shift  in  the  location  of  zero  pressure  as  discussed  above  has  also 
been  observed  experimentally58.  As  shown  by  comparing  f igures  6-1 2(a)  and  (b),  the 
effect  ot  the  foramen  magnum  opening  is  to  create  larger  tensile  pressures  for  longer  dura- 
tions at  anterior  points  in  the  brain  during  rear  impacts.  It,  therefore,  is  regarded  as  an  im 
portant  mechanism  in  head  injury  modeling. 


• I In  reverse  in  naturally  true  tor  the  coup  pressures,  the  coup 
pressure  s m the  open  skull. 


pressures  in  the  closed  skull  are  much  larger  than  the  coup 
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In  this  chapter  the  computed  dynamic  response  obtained  from  several  HIM  code  simula- 
tions m the  last  part  of  the  head  injury  model  development  program  are  presented.  They 
deal  with  the  human  configuration  and,  in  particular,  with  the  three-dimensional  discretiza- 
tion discussed  earlier  in  ( hapter  5 I hough  most  of  the  effort  described  throughout  this 
report  pertains  to  tin  development  of  a usable  tool  with  which  to  study  head  injury,  some 
of  the  results  presented  in  this  chapter  provide  useful  insight  into  the  mechanical  causes  of 
head  injure  an  important  objective  of  this  head  injure  model  study.  Fxccrpts  from  these 
results  were  presented  at  an  earlier  date  in  Reference  65. 

I rontal  Impact  Versus  Occipital  Impact 

I In  human  skull  model  employing  the  fine  discretization  (see  Figure  5 12)  was  sub- 
jected to  a typical  direct  impact  load  from  the  front  and  from  the  rear.  In  each  case  the 
peak  load  was  2,600  pounds  and  was  distributed  over  a skull  surface  area  of  about  two  sq  in. 
I he  impact  load  varied  with  time  as  a haversine  function  having  a duration  of  10  ms.  The 
direction  of  loading  nearly  passed  through  the  skull  center  of  gravity,  and  the  skull  was  com- 
pletely unrestrained,  fhe  skull  model  included  an  elastic  membrane  over  the  foramen  mag- 
num in  the  same  manner  as  previously  discussed  for  the  rhesus  skull  model.  Additionally, 
a compressible  bulk  modulus  of  305  psi  was  assigned  the  elements  constituting  the  subarach- 
noid space  I he  integration  time  step  size  was  1 ms. 

Cranial  Bone  Response  Minimum  and  maximum  principal  stresses  in  the  midsagittal 
plane  of  the  cranium  are  presented  in  contour  form  in  Figure  7-1  for  a frontal  impact.  Neg- 
ative and  positive  v alues  refer  to  compression  and  tension  stresses  in  the  skull  bone,  respec- 
tivclv  I or  the  frontal  load  condition,  the  outer  surface  of  the  frontal  bone  sustains  the 
maximum  compression  stress  (6,200  psi)  while  the  inner  surface  of  the  frontal  bone  sus- 
tains the  maximum  tension  stress  (2,050  psi).  These  are  bending  stresses.  They  appear  to 
accumulate  above  the  loading  region  where  the  frontal  bone  thickness  is  reduced.  The  rela- 
tive magnitudes  of  compression  and  tension  stresses  are  significant  since  tension  stress  is 
usually  thought  of  as  being  more  conducive  to  skull  fracture.  Yet  the  compressive  stresses 
are  roughly  three  times  greater.  If  they  were  equal,  then  fracture  would  be  expected  to 
begin  eventually  on  the  inner  and  outer  cranial  surface,  higher  up  on  the  frontal  bone.  But 
as  the  analysis  shows,  the  larger  compressive  stresses  may  significantly  alter  that  fracture 
pattern  by  initiating  fracture  on  the  outer  surface  of  the  frontal  bone  just  below  the  indi- 
cated position  for  impending  tensile  fracture. 

Principal  stress  data  for  the  occipital  impact  is  presented  in  Figure  7-2.  Maximum  com- 
pressive stress  (5,160  psi)  occurs  just  above  the  load  region  on  the  outer  skull  surface,  \gain 
these  are  bending  stresses.  The  maximum  tensile  stress  ( 1 ,820  psi),  however,  occurs  at  the 
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\ Iransient  structural  response  of  the  linear  skuN/hrain  system  in  Proceedings  ot  the  Nineteenth  Stapp 
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(a)  Compressive  stress. 
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I ipure  7 1 Principal  Stresses  in  Cranium  for  a I rontal  Impact  at  Peak  Response  Time. 


edge  of  the  foramen  magnum,  wet!  below  the  point  of  loading.  As  a result  the  foramen 
magnum  appears  to  present  an  area  of  stress  concentration  not  unlike  that  found  in  con- 
ventional engineering  shell  structures  in  the  vicinity  of  apertures.  I he  compressive  stresses 
are  again  almost  three  times  larger  than  the  tensile  stresses.  As  a result  their  contribution  to 
linear  fracture  cannot  be  ignored  even  though  tensile  stresses  may  be  more  conducive  to 
linear  fracture. 

The  effect  of  simulating  the  sandwiched  or  layered  nature  of  the  cranial  bone  structure 
is  illustrated  from  the  concentration  of  stress  contours  in  the  inner  and  outer  table  bone 
las  ers,  while  the  diploe,  or  middle  layer  of  softer  bone,  is  essentially  devoid  of  high  stresses. 

1 he  layered  construction  of  the  skull  is  indeed  very  efficient  in  resisting  the  bending  forces 
of  impact.  Compressive  bending  stresses  are  consistently  higher  than  tensile  bending  stresses 
in  the  above  two  simulations.  Because  bone  is  more  susceptible  to  failure  by  tension 
stresses,  it  is  difficult  to  pin-point  the  exact  local  of  impending  fracture.  Moreover,  no  ac- 
ceptable failure  criterion  exists  for  skull  bone  subjected  to  two-  and  three-dimensional  stress 
states.  It  is  clear  that  fracture  paths  and  areas  of  high  stress  are  considerably  altered  by  local 
geometrical  irregularities  of  the  cranial  structure.  It  is  also  clear  from  the  above  two  simula- 
tions that  the  occiput  dcvelopes  less  stress  for  the  same  load  than  does  the  frontal  bone  and 
apprears  to  be  stronger.  Linear  fractures  arc  apt  to  originate  on  the  inner  skull  surface  as 
well  as  the  outer  skull  surface,  f or  frontal  impacts  the  fracture  of  the  frontal  bone  is  likely 
to  occur  above  the  blow  at  a higher  position,  l or  occipital  impacts,  fracture  of  the  occipi- 
tal bone  is  likely  to  occur  either  above  or  below  the  impact  area  and  perhaps  as  far  down  as 
the  foramen  magnum. 

Intracranial  Response.  Data  in  this  section  as  in  all  HIM  code  predictions  are  presented 
for  the  midsagittal  plane  of  the  model  but  is  to  be  interpreted  as  predicted  data  for  a sagittal 
plane  just  to  the  right  or  left  of  the  midsagittal  plane  of  a human  skull.  I'his  is  because  the 
two  cerebral  hemispheres  are  actually  separated  by  the  falx  membrane  which  occupies  a 
good  portion  of  the  midsagittal  plane  A full  sagittal  section  of  cerebral  material  occurs  in 
actuality  only  to  the  right  or  left  of  this  plane. 

I he  results  for  each  load  case  arc  presented  together  in  figure  7-3  for  comparison.  Neg- 
ative and  positive  pressures  denote  compression  and  tension,  respectively.  Leak  pressure 
that  occurs  as  compression  in  the  anterior  brain  from  the  frontal  blow  (f  igure  7-3(a)  and 
(b))  is  less  likely  to  cause  brain  damage  than  the  peak  pressure  occurring  as  tension  in  the 
anterior  brain  as  induced  by  the  occipital  blow  (figure  7-3(c)  and  (d)).  They  arc  nearly  the 
same  in  absolute  magnitude,  but  the  brain  is  believed  more  susceptible  to  injury  from  ten- 
sion This  is  the  same  characteristic  a was  noted  with  the  rhesus  model  and  of  fers  a good 
explanation  for  the  statistically  observed  phenomenon  of  asymmetry  in  brain  injury  (see 
I able  6-4)  caused  by  frontal  and  occipital  impacts. 

Maximum  shear  strain  data  when  compared  for  the  two  load  cases  are  almost  indistin- 
guishable anil  therefore  appear  to  be  independent  of  whether  the  impact  is  frontal  or  occi- 
pital. I he  peak  shear  strain  for  the  occipital  impact  is  20%  greater  than  for  the  frontal  im- 
pact. Brain  damage  due  to  shear  strain  will  most  likclv  be  found  in  the  brainstem  or 
possibly  the  cerebellum  according  to  this  data,  regardless  of  the  impact  load  direction.  The 
midbrain  region  possesses  very  little  shear  strain. 
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Ihesc  intracranial  results,  both  pressure  and  maximum  shear  strain,  surely  indicate  that 
intracranial  pressure  would  be  the  mechanical  phenomenon  more  likely  to  explain  the  statis- 
tical distribuiton  ol  coup  and  contrecoup  damage  which  is  noted  in  Table  6-4.  I he  most 
severe  shear  strain  does  not  occur  at  either  a coup  or  contrecoup  location.  The  most  severe 
pressure  does  occur  at  these  positions.  The  location  of  the  most  severe  shear  strain  does  not 
v ary  with  the  load  direction,  but  the  location  of  the  most  severe  pressure  does. 

I Inis  tar  the  HIM  code  indicates  pressure  to  be  the  most  important  intracranial  response- 
tor  predicting  brain  damage.  However,  to  this  point,  only  translation  of  the  skull  has  been 
studied  and  whether  or  not  these  results  will  differ  for  rotation  remains  to  be  seen  The  ef- 
fect of  rotation  is  investigated  in  the  following  section.  ; 

Varying  the  Skull  Base  Boundary  Condition 

A series  ol  tour  HIM  code  simulations  are  presented  and  discussed  for  the  primary  pur- 
pose ot  determining  the  necessity  for  including  a neck  simulation  with  the  skull  model.  In 
tonducting  these  simulations  the  effect  of  skull  rotation  is  also  introduced.  Rigid  body 
motion  ot  the  skull  anti  the  response  of  the  brain  in  terms  ot  intracranial  pressure  data  and 
maximum  shear  strain  data  are  employed  in  the  discussion.  Cranial  bone  stresses  are  not 
discussed,  but  tor  each  simulation  like  those  of  the  previous  section,  information  pertaining 
to  linear  fracture  of  the  cranium  was  computed. 

I he  effect  ot  the  boundary  conditions  at  the  skull  base  upon  intracranial  response  is 
investigated  by  repeating  the  impact  simulation  for  a skull  which  is  (a)  fixed  against  move- 
ment at  the  base,  (b)  pinned  at  the  base  allowing  rotation  about  the  base,  (c)  rollered  at  the 
base  allowing  the  skull  to  slide  along  a plane,  and  (d)  completely  unrestrained  at  the  base- 
allowing  the  skull  to  translate  and  rotate.  I he  simulated  impact  load  was  a concentrated 
force,  varying  with  time  as  a haversine  function,  with  a duration  of  10  ms  and  a peak  force 
ot  1 ,000  pounds  | he  impact  load  is  applied  at  an  angle  of  1 5 degrees  with  the  horizontal 
on  the  frontal  bone  Again,  the  numerical  integration  time  step  size  was  1 ms. 

Skull  f ixed  at  the  Base.  The  fixed  skull  is  completely  restrained  at  each  of  two  node 
points  on  the  perimeter  ot  the  foramen  magnum  The  computed  displacements,  stresses, 
and  shear  strains  are  shown  in  figure  7-4  at  three  selected  times  corresponding  to  peak  load 
ing,  prior  to  peak  loading,  and  subsequent  to  peak  loading.  I he  sequence  of  displacements 
shows  very  little  rigid  body  motion  because  ot  the  restraint  at  the  base.  To  display  w hat 
rigid  body  motion  did  exist,  the  displacements  were  sealed  to  2 5 times  their  actual  values. 

In  the  subarachnoid  region  near  the  top  of  the  skull  a relative  motion  between  the  inner 
skull  elements  and  the  outer  brain  elements  can  be  observed  at  5 ms  (note  the  |og  at  the  top 
along  the  vertical  mesh  lines).  I lus  is  graphical  evidence  of  the  desired  effect  designed  by 
specifying  a relatively  compressible  bulk  modulus  in  the  subarachnoid  elements  I his  result 
can  be  observed  in  other  HIM  code  simulations  to  a lesser  or  gn. .iter  extend  I hough  the 
relative  motion  is  always  slight,  it  serves  to  reduce  significantly  the  stresses  and  strains  in 
the  brain  that  would  otherw  ise  develop  by  assuming  that  the  brain  is  riguilv  f ixed  to  the 
inner  skull  surface. 

I he  curvature  ot  the  stress  contours,  for  which  values  are  given  in  I able  7-1 . demon 
strates  the  rotation  taking  place  to  the  virtual  exclusion  of  translation  \ significant  amount 
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Figure  7-4  Direct  Impact  of  l ixcd  Skull. 


ol  tension  is  evident  in  the  brai  < in  tins  simulation.  I'he  skuli/is  aceelerateil  “hack  over"  the 
support,  in  it  tails  causing  compressi  m in  the  anterior  and  tension  in  the  posterior  brain  I he 
skull,  however,  soon  accelerates  back  in  the  opposite  direction  as  the  load  subsides  I bis 
causes  a sign  reversal  in  the  stresses  of  the  anterior  brain.  Actually  the  peak  tension  stress 
occurred  at  9 ms  (not  shown)  on  the  anterior  brain  surface  in  this  simulation.  I'he  elasticity 
ol  the  skull  bone  has  participated  significantly  because  of  the  restraining  condition. 

Shear  strains  develop  initially  in  the  brainstem  region  then  subside  at  7 ms,  only  to  in- 
crease again  to  the  maximum  v alue  at  9 ms.  I bis  oscillation  is  a manifestation  of  the  same 
phenomenon  discussed  above  for  the  stress  response.  I here  is  only  a slight  tendency  for  the 
shear  strains  to  develop  on  the  outer  surface  of  the  brain.  I he  brainstem  sustained  virtually 
all  the  shear  strain.  It  might  reasonably  be  expected  that  this  is  due  to  the  fixed  support 
condition  but  results  will  be  shown  for  an  unrestrained  skull  which  discount  that  possibility 

Skull  Hinged  at  the  base.  Only  in  the  most  rare  situation  would  a completely  restrained 
skull  simulation  be  appropriate.  A more  realistic  constraint  would  allow  some  rigid  body 
movement  of  the  skull  The  results  for  a hinged  skull  are  shown  in  f igure  7-5.  I he  rotation 
takes  place  about  a pinned  node  on  the  anterior  ridge  of  the  foramen  magnum  and  can  be 
termed  noncentroidal  rotation.  Significant  rigid  body  motion  occurs  as  compared  with  the 
completely  restrained  skull  and  as  a result  the  intracranial  response  is  very  different,  f irst 
of  all  there  is  no  reversal  of  sign  in  the  computed  intracranial  response  and  second,  the  com- 
pressive stresses  and  shear  strain  are  considerably  larger  as  shown  by  the  values  in  fable  7-2 
Ibis  condition  is  therefore  more  severe  ’ban  the  completely  restrained  condition  with  regard 
to  intracranial  response.  1 his  is  not  true  for  the  skuII  bone  response. 

I he  curvature  ol  the  stress  contours,  especially  in  the  anterior  brain,  reflect  the  rotation 
(rigid  body)  that  is  taking  place.  I he  contrecoup  response  is  comparatively  mild  as  ev  idence 
in  the  peak  response  at  5 ms  showing  less  than  1 3-ps>  tension  in  the  posterior  brain. 

The  maximum  shear  strain  response  peaks  at  7 ms;  i.».\,  subsequent  to  peak  pressure  re- 
sponse The  brainstem  region  sustains  the  more  severe  shear  strain  and  the  midbrain  region 
sustains  the  least  shear  strain  Although  not  shown,  these  observations  are  also  true  of  the 
computed  normal  strain  distribution. 

Skull  Sliding  on  Its  Base.  In  this  instance  the  simulated  skull  is  allowed  to  slide  on  a 
horizontal  plane  tangent  to  the  skull  base  under  the  influence  of  the  impact  force  I he  rigid 
body  motion  in  this  case  is  primarily  translational. 

t he  intracranial  response  is  basically  similar  to  the  hinged  skull  condition  as  shown  in 
figure  7-6  I he  stress  contour  values  presented  in  fable  7-3  indicate  a higher  compressive 
stress  in  the  anterior  brain  and  also  a slightly  more  severe  contrecoup  condition  in  the  pos 
terior  brain  I he  magnitude  of  tensile  stresses  are  higher  ami  distributed  over  about  the 
same  amount  of  brain  material  as  in  the  hinged  skull. 

I'he  shear  strain  contour  response,  with  regard  to  distribution,  looks  tin  same  as  the  re- 
sponse for  the  hinged  skull  However  the  magnitude  of  these  contours  is  lower  for  the  slid 
mg  skull  that  it  is  for  the  hinged  skull.  So,  while  the  model  predicts  the  translated  skull  to 
be  more  severe  with  regard  to  coup  and  contrecoup  stresses,  it  also  predicts  the  rotated  skull 
to  be  more  severe  regarding  the  strain  condition  in  the  brainstem  region 


I 33 


Tabic  7-2.  Contour  Values  for  Hinged  Skull 


( iontour 
No. 

Stress3 

(psi) 

Maximum 
Shear  Strain 

1 

-40.8 

0.093 

2 

34.0 

0.184 

3 

-27.1 

0.276 

4 

-20.3 

0.367 

5 

13.5 

0.459 

6 

-6.7 

0.550 

7 

+0.2 

0.641 

8 

+ 7.0 

0.733 

9 

+ 13.8 

0.824 

10 

+20.6 

0.916 

1 Negative  sign  indicates  compression,  and  positive  sign 
indicates  tension. 
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I aide  7-3  Contour  Values  for  Sliding  Skull 


a Negative  sign  indicates  compression,  and  positive  sign 
indicates  tension. 
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I igurc  7 6 Direct  Impact  of  .1  Muling  Skull 


Skull  Completely  I •"ret*  at  the  Base.  I he  response  ol  the  brain  is  slightly  more  severe 
than  previous  responses  it',  tor  the  same  impact  lo.ul.  the  boundary  restraints  at  the  base  of 
the  skull  are  removed.  I he  HIM  eode  results  lor  the  impacted,  free  skull  arc  presented  in 
Figure  7-7  Here  rigid  body  motion  is  significant  and  is  a general  or  curvilinear  state  of  mo- 
tion combining  both  rotation  and  translation.  Response  parameters  are  presented  at  the 
same  selected  times  as  before. 

The  stress  contour  lines  are  generally  more  straight  and  less  vertical  than  in  the  previous 
eases.  I his  does  not  necessarily  mean  less  rotation  occurred,  but  only  that  there  is  in  this 
case,  more  of  an  influence  from  translation.  I he  rotational  influence  is  manifested  in  the 
change  in  orientation  of  the  rather  straight  contours.  I he  stress  magnitudes  reported  in 
Table  7-4  show  that  higher  contrecoup  brain  stresses  occurred  in  the  free  skull  than  in  any 
of  the  previous  simulations.  I he  stress-time  response  follows  closely  the  load-time  specifi- 
cation with  the  peak  response  occurring  at  5 ms.  I he  influence  of  the  foramen  magnum  is 
evident  from  the  posterior  location  of  the  zero  pressure  region.  These  two  observations  are 
generally  true  of  all  these  simulations. 

Sheat  strains  develop  initially  in  the  brainstem  region  as  before.  They  continue  to  in- 
crease throughout  the  simulation  to  their  maximum  value.  I he  maximum  response  time  is 
9 ms  after  which  they  subside.  At  any  given  time  the  shear  strain  diminishes  with  distance 
from  the  foramen  magnum  towards  the  midbrain  region  where  the  shear  strain  prediction  is 
always  observed  to  be  minimal.  Beginning  from  the  midbrain  region,  the  shear  strain  gener- 
ally increases  with  radial  distance  toward  the  outer  surface  of  the  brain,  particularly  the 
superior  surface.  The  shear  effect  can  also  be  seen  from  the  displacements  in  the  distortion 
of  the  outer  brain  elements  as  compared  with  the  distortion  of  the  inner  brain  elements. 
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I able  7 4 C ontour  Values  for  f ree  Skull 


ontour 

No. 

Stress2, 

(psi) 

Maximum 
Shear  Strain 

i 

-51.9 

0.091 

2 

-42.0 

0.180 

3 

-32.1 

0.271 

4 

-22.2 

0.361 

5 

-12.3 

0.451 

6 

-2.46 

0.541 

7 

+7.43 

0.631 

8 

+ 17.3 

0.721 

9 

+27.2 

0.812 

10 

+ 37.1 

0.902 

^Negative  sign  indicates  compression,  and  positive  sign 
indicates  tension. 

**  1 psi  = 0.069  bar. 


8 INDIRl  C I IMP  AC  I AC<  I I I K A I ION 


In  many  hazardous  dynamic  environments  occupants  are  restrained  so  that  their  heads 
may  not  directly  impact  surrounding  objects.  Yet  there  is  a need  to  determine  the  mechan- 
ical stress  and  strain  response  of  the  brain  under  these  circumstances.  I h is  is  also  true  of 
human  volunteer  testing  where  the  environment  should  not  be  truly  hazardous,  but  mech- 
anical stress  and  strain  predictions  are  necessary  for  correlation  with  measured  physiologi- 
cal response  In  these  instances  the  impact  forces  arc  transmitted  to  the  head  at  the  base  of 
the  skull  and  subject  the  head  to  impact  acceleration  as  described  by  Kwing  and  Thomas06. 

An  intensive  effort  is  ongoing  to  quantify  the  human  dynamic  response  to  impact  accel- 
eration and  its  correlation  with  injury  and  physiological  effects.  This  research  is  being  con- 
ducted at  thi  Naval  \crospacc  Medical  Research  Laboratory  (NAMKL),  Michoud  Station, 
New  Orleans,  Louisiana.  Volunteer  subjects  and  human  analogs  are  given  impact  accelera- 
tions at  a modern  pneumatically  powered  sled  facility.  A computer-controlled  data  acqui- 
sition system  continuously  monitors  all  input  parameters  as  well  as  selected  output  param- 
eters of  the  subject's  response.  Acquired  data  is  then  applied  in  the  development  of  mathe- 
matical models  which  are  intended  to  complement  and  extend  knowledge  acquired  experi- 
mentally Once  validation  of  these  models  is  achieved  and  once  sufficient  level  of  confi- 
dence has  been  attained  with  the  models,  they  will  play  an  active  role  in  the  design  and 
evaluation  of  protective  systems  for  humans  in  hazardous  dynamic  environments.  Specifi- 
cally, alternative  aircraft  anil  ship  designs  can  be  evaluated  from  the  standpoint  of  crew 
system  safety,  with  the  mathematical  models  prov  iding  basic  data  in  the  evaluation  process. 
It  is  never  too  soon,  when  considering  the  development  of  new  high  performance  craft,  to 
ask  whether  crews  can  function  during  normal  operation  and  be  reasonably  protected  during 
abnormal  operation  of  high  performance  craft.  Hydrofoils  and  surface  effect  machines  are 
intended  to  travel  at  speeds  up  to  80  knots  at  rough  sea  states. 

I o establish  the  HIM  code  as  a usable  tool  in  instances  of  indirect  impact,  several  pre- 
liminary problems  required  solutions.  ( I he  necessary  theoretical  considerations  were  pre- 
sented in  ( ihapter  3. ) 

Mass  Distribution  Properties  of  the  Skull 

In  previous  studies  with  the  finite  element  head  injury  model,  it  was  found  that  signifi- 
cant intracranial  pressures  occur  well  after  skull  deformation  modes  had  dissipated.  I his 
suggests  the  importance  of  computing  the  correct  rigid  body  response  of  the  head  due  to 
dynamic  loading.  It  is  therefore  necessary  that  the  discretized  finite  element  model  of  the 
head  possess  the  correct  inertial  properties.  I o this  end  the  capability  for  computing  the 
mass  moment  ol  inertia  tensor  w ith  respect  to  the  center  of  gravity  was  added  to  the  exist 
ing  HIM  code  preprocessor  I ventually  this  computed  data  was  to  be  compared  with  mea- 
sured inertial  distribution  properties  of  rhesus  skulls  used  in  NAMRI.  testing.  Correlation 
would  validate  the  accuracy  of  the  discretization. 


Naval  V rnspui  c Me.  I it  .(I  Research  laboratory  N VMKI 
jiirkr.it mti  by  ( harming  I I sung  and  Manirl  | I bonus 


Monograph  ’I  Human  head  and  ruck  response  to  impact 
I'cnsacola  I I , Aug  I ‘>7? 
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Measurement  of  these  properties  is  difficult  but  has  been  achieved  for  human  cadaveric 
skulls'’  Numerical  computation  of  the  inertia  tensor  for  the  modeled  skull  is  conceptually 
very  easy,  however,  for  example,  the  mass  moment  of  inertia  tensor  component  with  re- 
spect to  the  \ -axis  ly  through  the  skull  is  defined  as 


/<*2 


/^Idm 


where  the  integral  is  evaluated  over  the  entire  skull.  Since,  in  the  finite  element  discretiza- 
tion process,  the  skull  continuum  has  been  conveniently  divided  into  a great  many  small 
elements,  I . can  be  approximated  as 

e 


where  xj  and  z.j  are  the  distances  in  the  x and  z directions  from  the  y-axis  to  the  centroid  of 
the  ith  element,  mj  is  the  mass  of  ith  element  and  e is  the  number  of  elements  constituting 
the  discretization.  The  element  mass  mj  is  computed  numerically  by  integrating  the  element 
volume  which  in  turn  is  proportional  to  the  determinant  of  the  Jacobian  j J | , which  was 
defined  in  the  element  stiffness  relationship  (section  3.0).  1 his  quantity  is  routinely  evalu- 
ated for  each  clement  during  the  initial  check  of  the  finite  element  mesh  subsequent  to  the 
automatic  mesh  generative  process  and  prior  to  HIM  code  execution.  So  the  most  impor- 
tant information  for  computing  inertial  quantities  is  available  and  little  additional  effort  is 
required.  It  is  also  necessary  to  code  the  parallel  axis  theorem  to  enable  evaluation  of  the 
inertia  tensor  with  respect  to  various  head  anatomical  reference  systems. 

Results  of  the  mass  moment  of  inertia  computations  are  presented  in  fable  8-1  for  the 
rhesus  skull  model  used  in  the  validation  efforts  reported  in  Chapter  6.  I he  computational 
scheme  was  checked  by  measuring  the  displaced  volume  of  an  empty,  dry  rhesus  skull  and 
successfully  comparing  it  with  the  computed  volume  of  a discretized  version  of  the  same 
skull  In  the  absence  of  measured  inertial  quantities,  no  direct  check  of  the  data  in  fable 
8-1  is  possible.  But  the  results  appear  reasonable  and  also  compare  well  with  the  inertia  of 
a solid  sphere  about  a diametrical  axis  that  is  equivalent  to  the  rhesus  skull  in  total  mass. 

I he  computational  method  does  depend  on  having  a great  many  small  elements  for 
accuracy. 


Vtv.il  Aerospace  Medical  Research  l aboratory.  VVMKI  1193  Measurement  of  mass  distribution  parameters  of  ana 
iiimiv.il  segments,  by  I dvsard  B Becker  fensacola,  I I , Oct  1973. 


142 


frequence  Response 


Frequency  response  data  ot  r lie  head  injury  model  is  usef  ul  for  two  different  reasons 
and  such  data  should  be  obtained  in  connection  with  validation.  I list,  the  data  will  aid  in 
the  selection  of  appropriate  frequency  response  ranges  for  transducers  used  in  experimental 
validation  tests  ft  w ill  provide  upper  and  lower  bounds  to  the  spectrum  of  dominant  fre 
quench  , possessed  by  the  head. 

No  less  important  is  the  utility  of  frequency  response  data  regarding  selection  of  appro 
pi  'ate  time  step  size  in  numerical  simulations.  Because  the  effect  of  discretizing  the  time 
domain  during  any  execution  of  the  model  is  to  numerically  filter  the  "true”  solution,  it  is 
important  that  the  ''residual”  or  computed  solution  include  the  major  frequency  compo- 
nents fo  achieve  this  condition,  the  selection  of  a time  step  size  when  executing  the  HIM 
code  must  be  based,  at  least  in  part,  on  f requency  response  information. 

I requenev  components  tor  various  1 1 1 VI  code  models  were  determined  by  using  a gen 
cral  structural  analysis  computer  code  called  SAP,,X  . I his  code  extracts  the  eigenvalues 
from  large  structural  stiffness  matrices  bv  a technique  known  as  subspace  iteration6'* . I lie 
method  is  the  current  recommended  procedure  when  large  matrices  arc  involved.  Yet,  there 
is  some  doubt  as  to  whether  subspace  iteration  converges  to  the  lowest,  and  most  important, 
mode.  In  any  ease  the  user  of  this  method  must  accept  the  numerical  results  with  little  as- 
surance beyond  a convergence  to  a prescribed  tolerance  that  the  lowest  eigenvalue  is  realis- 
tically the  fundamental  frequency, 

1 he  natural  frequencies  associated  with  the  first  four  modes  of  the  HIM  code  rhesus  and 
human  skull  models  are  presented  in  labie  8-2.  I Experimental  measurements  of  the  lowest 
natural  frequency  has  yielded  numbers  which  arc  in  the  range  of  1 to  20  Hz.  for  human  and 
rhesus  skulls.  So  the  computed  values  do  appear  high  I his  may  be  attributed  to  the 
model's  skull  bone  stiffness  which,  from  static  validation  efforts  (Chapter  6),  was  shown  to 

be  toi,  high. 


( onversion  of  N AMR  I Measured  Data 

A large  data  base  constructed  and  stored  at  NAMR1.  contains  the  measured  response  of 
the  head  and  neck  for  a great  many  dynamic  tests  of  humans  and  human  analogues  I he 
crucial  data  conveying  the  motion  of  the  head  was  extracted  from  this  source  and  applied  as 
input  to  the  finite  element  head  injury  model.  In  this  manner  the  intracranial  response  pre 
dieted  by  the  model  would  be  consistent  with  input  motion  ami  would  also  include  the  in 
fluenee  of  the  neck.  A key  issue  is  the  technique  in  which  the  measured  motion  is  applied 
over  the  eontiuum  of  tin-  head,  that  is,  the  means  by  which  the  motion  is  distributed  to 
drive  the  skull  model. 


I m versify  of  ( alitnrma.  i iillcgr  of  I nginrering  Kt-port  I I l<(  7J  I I S \|»  IV,  a struniir.il  analysis  program,  for  stain 
and  dynamic  response  of  linear  systems.  b\  K | Bathe,  I I Wilson,  and  I I Peterson  Berkelm  . < V Jun  197$ 

69 

K | Batric  and  I I Wilson  Solution  methods  for  eigenvalue  problems  m structural  mechanics.  International 
journal  of  \ micrual  Methods  m I ngincmng.  vol  6.  1973.  pp  - 1 3 22<* 
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I' able  8-2  ( omputed  Natural  Frequencies  lor  HIM  (.ode  Models 


1 

Mode 

frequencies  (II/.) 

1 mpty 
Rhesus 

Full 

Rhesus 

Human 

~ 

1 

130 

108 

43 

2 

838 

120 

51 

3 

1.121 

126 

55 

4 



2,266 

127 

57 

I he  theory  in  ( hapter  3 demonstrates  the  feasibility  of  prescribing  displaeenu  nt  histor- 
ies of  selected  node  points  on  the  skull  model.  I he  method  ol  prescribing  displ.u  ement  his 
tories  was  selected  at  the  outset  because  the  current  version  ol  the  HIM  code  did  nm  . , : 1 < « >a 
lor  specification  of  acceleration  histories  nut  purported  to  work  with  prescribed  displ.n 
merit  histories  and  initial  velocities  as  well. 

I he  following  describes  how  the  measured  displacement  data  is  transformed  to  provide 
prescribed  input  displacement  history  data  for  the  HIM  code  model'"  VYMRI  data  which 
is  supplied  on  magnetic  tape,  includes  the  translation  ol  the  origin  ol  the  anatomical  coor 
dilute  sv  stem  for  the  head  as  the  subject  is  carried  along  on  the  sled  \r  selected  sample 
times,  the  three-dimensional  rotation  of  the  head  anatomical  coordinate  system  is  supplied 
along  with  the  translation  of  the  system.  I he  rotation  data  consists  cither  of  tour  qu.iter 
mons  or  three  luilcr  angles  I hus,  the  displacement  history  of  each  node  point  on  tin  skull 
brain  model  selected  for  driving  is  obtained  by  the  following  translorm.it i<>: 


Vervonal  communication  with  I S I untie . Naval  Aerospace  Kent  arch  I al>ora(or\  Detachment  Wtchoutl  Station  New 
Orleans.  I A.  I ft  Sep  1975 
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when-  the  subscript  zero  denotes  the  location  of  the  anatomical  system  origin  in  laboratory 
coordinates  at  time  t and  the  subscript  s denotes  the  location  of  a selected  node  point  on 
the  model  in  the  anatomical  coordinate  system.  I he  term  | I (t)  | is  the  transformation 
matrix  from  the  anatomical  coordinate  system  to  the  laboratory  coordinate  system  at  time  t 
and  is  defined  in  I able  8-3 . 

Because  the  HIM  code  does  not  allow  a variable  time  step  specification,  the  transformed 
displacement  histories  must  be  for  a constant  time  interval.  I his  presents  no  problem  when 
sensor  data  is  transformed  because  it  is  currently  supplied  at  equally  spaced,  1-ms  intervals. 
However,  il  the  data  is  supplied  from  high  speed  photography  it  is  currently  only  approxi- 
mately cquallv  spaced  (2  ms  intervals)  and  interpolation  will  then  be  neccssarv  Since  f iller 
angles  can  be  discontinuous  they  are  unsuited  to  interpolation;  quaternion  would  be  recom- 
mended in  cases  where  interpolation  is  necessary. 

Because  \ WIKI  data  acquisition  is  completely  three-dimensional,  it  senses  and  records 
the  inevitable  movement  ol  the  skull,  which  is  not  symmetrical  with  the  midsagittal  plane 
1 he  HIM  code  can  provide  a completely  three-dimensional  discretization,  but  a half-skull 
discretization  is  being  employed  in  this  effort  for  economical  reasons.  I bis  means  that  any 
data  contributing  to  unsymmctrical  displacement  must  be  filtered  in  the  transformation 
process  I his  is  relatively  easily  accomplished  by  zeroing  out  the  first  and  third  I ulcr  angles 
prior  to  entering  into  the  transformation  equation,  along  with  the  v translational  displace- 
ment component  after  going  through  the  computation.  I he  computed  angular  displacement 
of  the  head  in  the  midsagittal  plane  is  shown  in  f igure  8-1 

1 he  selection  ol  node  points  whose  displacement  histories  are  to  be  prescribed  requires 
some  elaboration.  Basically,  a minimum  ol  three  noncollinear  nodes  should  be  selected,  and 
they  should  be  on  the  skull  bone  Since  the  transformation  computation  mav  involve  some- 
error,  no  matter  how  slight,  the  prescribed  displacements  for  several  nodes  will  not  be  exact- 
ly consistent  w ith  rigid  body  motion.  I hercforc.  a fictitious  strain  will  have  been  intro- 
duced into  the  skull  bone  and  will  depend  on  the  error  magnitude.  I his  difficulty  can  be 
minimized  In  selecting  a minimum  number  of  nodes  and  nodes  which  are  physically'-  located 
at  large  distances  from  one  another  and  also  by  carrying  out  the  prescribed  displacement 
- imputation  to  the  maximum  number  ol  decimal  places  allowed  by  computer  word  length 
I he  extent  to  which  the  fictitious  skull  strain  affects  the  brain  response  is  believed  to  be 
negligible  And  since  the  brain’s  response  in  indirect  impact  simulations  is  w hat  is  primarily 
sought,  the  significance  of  the  fictitious  skull  strain  is  mitigated. 

Prcliminarv  Simulations 

\ m liminary  stud x of  t he  HIM  code's  ability  to  operate  with  prescribed  displacements 
w i a luted  with  a simple  model  (I  igurc  8 2),  but  one  which  retained  the  basic  ingredients 
ol  i l ull  model  I his  study  is  important  in  revealing  the  differences  in  numerical  behav  ior 
between  solutions  for  direct  and  indirect  impact. 

I irst  tin  model  was  sub|cctcd  to  a prescribed  forcing  function  (direct  impact),  and  the 
resulting  displacements  were  plotted.  A curve  fit  to  this  displacement  data  yielded  a second- 
order  poly  nomial  displacement  history  as  shown  in  f igure  8 3(a).  Next,  this  displacement 
history  was  prescribed  in  a second  execution  of  the  simple  model  (with  a zero  forcing  func- 
tion this  time)  I he  “intracranial  pressure"  results  from  both  runs  compared  reasonably 


Table  8-3.  Definition  of  Transformation  Matrix  Coefficients  in  Terms  of  Euler  Angles  or  Quaternions 
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Pressure  (psi) 

Pressure  in  First  Fiement  <psi> 

fension  Compression  tension  Compression 


I 


well  as  is  seen  in  figure  8 3(b)  anil  (e)  considering  the  differences  in  input.  I his  indicated 
that  the  specified  displacement  option  of  the  HIM  code  was  apparently  working  satisfac- 
torily. 

f ollowing  this,  further  preliminary  computations  were  conducted  with  two  oscillating- 
type  functions  (sine  squared  and  sine)  having  longer  durations.  I his  type  of  displacement 
function  was  anticipated  with  NAMRI.  data  extracted  from  both  repeated  load  (obtained 
from  shake  table  experiments)  and  sled  tests,  flic  input  function  and  the  results  are  pre- 
sented in  f igures  84(a)  and  (b).  Computed  “intracranial  pressures”  in  each  case  fluctuated 
unexpectedly  Suspecting  that  the  time  step  size  was  to  coarse  in  comparison  to  the  natural 
period  of  a typical  brain  element,  the  time  step  was  successively  reduced  first  for  the  sine- 
squared  function  and  then  for  the  sine  function.  In  f igure  84(a)  it  can  lie  seen  that  the 
time  step  size  reduction  had  no  effect. 

I he  initial  conditions  were  then  investigated.  I'he  sine  squared  function  has  a zero  first 
derivative  at  time  zero  (initial  velocity ),  but  a nonzero  second  derivative  at  time  zero  (initial 
acceleration).  By  assigning  these  initial  conditions  and  repeating  the  computer  run  with  the 
most  coarse  time  step,  the  computed  “intracranial  pressures”  were  found  to  varv  smoothly 
in  an  expected  manner.  However,  a paradox  arises  and  can  be  seen  by  writing  the  equation 
of  motion  for  a simple,  undamped,  spring-mass  system  (see  liquation  3-30).  1'his  equation 
is  obviously-  not  satisfied  at  time  zero  when  an  initial  acceleration  and  a zero  initial  displace- 
ment is  applied,  finis,  while  the  sine-squared  function  for  displacement  definitely  possesses 
a nonzero  second  derivative  at  time  zero  and  through  its  specification  the  results  appear 
satisfactory,  it  is  not  strictly  valid  to  specify  it  as  an  initial  condition  for  the  physical 
system. 

I he  sine  function  for  displacement  possesses  a nonzero  first  derivative  (initial  velocity) 
and  a zero  second  derivative.  Specification  of  these  initial  conditions  reduces  considerably 
the  computed  fluctuation  in  pressure  as  is  seen  in  f igure  84(b).  f urther  reductions  in  time 
step  size  then  show  still  more  improvement,  but  this  improvement  would  not  exist  without 
first  having  specified  the  correct  initial  velocity. 

Therefore,  it  is  necessary  to  be  cognizant  of  initial  velocity  conditions  and  specify  them 
when  using  the  displacement  option  of  the  HIM  code.  It  is  also  apparently  necessary  to 
employ  time  steps  no  larger  than  eight  times  the  critical  value  to  avoid  excessive  intracranial 
pressure  fluctuations.  Since  ty  pical  NAMRI.  sled  tests  have  durations  of  around  300  ms,  the 
cost  of  using  the  present  operator  will  probably  be  high.  It  may  be  possible  to  use  the  same- 
operator  (New  mark  f>  method)  along  with  some  prescribed  numerical  damping.  The  fluctu- 
ations which  appear  in  the  pressure  solution  are  not  unstable  (they  do  not  grow  without 
bound).  It  is  seen  that  their  mean  or  average  value  is  the  solution  sought,  and  thus  some 
form  of  damping  flj  4 1/4)  could  allow  the  use  of  larger,  more  economical  time  steps. 

HIM  Code  Simulation  of  Indirect  Impact 

Using  the  human  skull  model,  several  attempts  were  made  to  compute  intracranial  pres- 
sure histories  while  specifying  input  displacements.  In  each  case  the  intracranial  pressure 
results  exhibited  the  fluctuation  noted  in  the  previous  section.  Durations  of  these  pre- 
scribed displacement  histories  varied  from  10  to  100  ms  and  time  step  sizes  were  from  1 t 
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5 ms  A tvpic.il  example  is  illustrated  m l igurc  8 5 Here  tlu  .nput  disp>.w<  m<  tit  1 1 iguu 

8 5a)  is  translational  onK  and  in  tins  ease  is  t lie  measure  I d ispi  it  i n • • 1 1 ' e m> 

the  initial  phase  of  NAM K I sled  test  no  1X0567  I hese  displa.  'its  u<  rt 

assigned  to  several  node  points  surrounding  the  foramen  magnum  s. . tun  tin  si  nil  m.idci  w is 
driven  at  the  base  of  the  skull.  I'ranslation  occurred  in  the  < . md  (,  dn<  i ta>m  simul 
taneouslv . As  can  he  seen  in  Figure  8 1 , little  rotation  of  the  head  on  urud  during  the  lust 
100  ms  of  run  no.  I.X0567,  fhc  pressure  fluctuations  taken  from  an  eh  ment  n tin  jmsiu 
u>r  brain  are  shown  in  Figure  8-5(l>).  I his  behavior  was  exhibited  throughout  tin  hi  un 
I he  t ime  step  si/e,  5 ms,  is  believed  to  he  too  coarse  lor  use  with  the  brain  eh  rents  w hu h 
are  similar  in  size  to  those  for  the  preliminary  model  diseussed  in  the  previous  v .non  N et 
a smaller  time  step  would  not  he  economical  lor  a NAMRI.  test  duration  ol  3o<>  r I lus 
anah  sis  cost  about  $200  and  spanned  only  the  lirst  one-third  of  the  total  test  i un.  I has.  a 
more  economical  integration  scheme  is  required  that  will  a I loss  larger  time  steps  and  contain 
sufficient  numerical  damping  to  produce  smooth  computed  pressure  histories  m the  brain 
Ihe  smoothed  average  of  the  pressure  fluctuation  exhibits  a tensile  pressure  in  the  posterior 
brain  as  would  be  expected  f urther,  the  pressure  magnitude  is  considerable  lower  than 
HIM  code  pressure  predictions  lor  direct  impact  (Chapter  7).  I his  appears  to  he  uniformly 
true  with  comparisons  of  indirect  and  direct  impact  intracranial  response 

A possible  alternative  to  specifying  directly  an  input  displacement  is  to  compute  before 
hand  an  equivalent  forcing  function  which  when  specified  would  produce  the  desired  dis- 
placement I he  computed  inertial  properties  (Table  8-1  ) of  the  model  must  be  employed 
m this  process  to  transform  a prescribed  displacement  to  an  "equivalent'’  prescribed  forcing 
function. 

I o demonstrate  this  method  the  nodal  points  on  the  perimeter  of  the  foramen  magnum 
were  assigned  equivalent  force-time  histories  and  a simulation  was  made  I hese  nodes  were 
constrained  to  translate  in  the  (i  direction  according  to  the  displacement  function  U(t) 

500  t + 1200  t , where  U is  measured  in  inches  anti  t in  seconds,  I hcrcforc,  the  skull  moves 
according  to  an  initial  velocity  of  500  in./s  and  a constant  acceleration  of  about  6.2  g's. 

I very  node  in  the  skull-brain  system  was  automatically  assigned  the  initial  velocity  value. 
Initial  displacements  were  assumed  zero,  although  provision  is  made  lor  their  simultaneous 
specification.  The  simulation  was  carried  out  through  10  ms 

The  displacement,  stress,  and  shear  strain  are  presented  in  figure  8-6  at  5,  7,  and  9 ms. 

I he  stress  state  continues  to  increase  until  at  9 ms  the  maximum  stresses  are  developed.  A 
slight  forward  rotation  is  evidenced  in  the  curvature  of  the  stress  contours  of  the  anterior 
brain.  It  cannot  be  seen  in  the  displacements  because  the  translational  movement  is  dome 
mint.  I he  stress  and  strain  magnitudes  are  given  in  I able  84  and  show  lower  magnitudes 
than  those  in  the  direct  impacts,  which,  of  course,  involve  higher  acceleration  levels 

Linear  .Model  Limitation 

It  should  be  emphasized  that  the  HIM  code  is  a linear  finite  element  code  and  that  all 
the  results  discussed  to  this  point  are  subject  to  the  limitations  imposed  by  linear  elasticity 
theory  I wo  such  limitations  which  are  particularly  important  involve  geometry  and  are  the 
limitations  of  small  displacements  and  small  strains.  The  normal  strain  component  i x ol 
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l igure  8-6  Indirect  Impact  Response  tor  translated  Skull 


I able  8-4.  Contour  Values  for  Translated  Skull 


Contour 

No. 

Stress1, ^ 
(psi) 

Maximum 
Shear  Strain 

1 

-9.0 

0.029 

2 

-7.0 

0.059 

3 

-5.0 

0.088 

4 

3.0 

0. 1 1 7 

5 

4). 9 

0.147 

6 

+ 1.1 

0.176 

7 

+ 3.1 

0.205 

8 

+ 5.2 

0.233 

9 

+ 7.2 

0.263 

10 

+9.2 

0.293 

•'Negative  sign  indicates  compression,  and  positive  sign 
indicates  tension. 

1 psi  0.069  bar. 


I (.{nation  3 4 is  defined  .is  t.w  change  in  length  ol  a unit  length  originally  oriented  parallel 
to  the  X axis  I he  shear  strain  component  c xy  of  liquation  3-5  is  defined  as  half  the  angle- 
change  between  two  lines  originally  parallel  to  the  X and  V axes  I he  second -order  strain 
terms  ot  1, (.{nations  3-4  and  3-5  extend  these  delinitions  to  cases  involving  large  movements 
which  have  rotated  the  original  axes  bv  gross  amounts. 

\ stud v w as  made  to  determine  the  extent  of  rotation  that  the  HIM  code  simulation 
could  sustain  before  the  linear  assumption  was  no  longer  valid.  I vidence  of  error  becomes 
apparent  when  the  computed  displacements  exhibit  a uniform  dilation.  The  studv  consisted 
merely  ol  trving  several  head  impact  simulations  allowing  rotation  to  occur  and  measuring 
the  angle  of  rotation  ( icncrally  it  was  found  that  dilation  in  the  displacements  began  occur- 
ring at  about  10  degrees  of  rigid  body  rotation.  The  dilation  is  uniform  and  is  not  accom- 
panied In  the  development  of  any  stress  or  strain  In  fact,  the  stress  and  strain  responses 
appear  to  be  normal,  i c.,  the  response  as  usual  subsides  subsequent  to  application  of  peak 
loading  even  though  the  rigid  rotations  continue  beyond  10  degrees  Still  the  validity  of  the 
linear  computation  must  be  regarded  as  suspect  for  large  rotations. 

I o correct  the  situation,  a nonlinear  analysis  is  required.  One  method  is  to  update  the 
coordinate  system  for  the  skull-brain  system  during  the  latter  time  steps  when  large  rota- 
tions are  being  experienced.  I his  is  theoretically  an  easy  modification  in  the  HIM  code  and 
does  not  require  the  addition  ol  the  second-order  terms  in  the  strain  formulation.  After 
each  time  step  the  computed  displacements  are  added  to  the  coordinates  to  define  a new 
system  Naturally  it  is  costly  computationally,  latch  time  that  the  coordinate  system  is  up- 
dated the  modified  stiffness  matrix  on  the  left-hand  side  of  liquation  3-34  will  change,  nec- 
essitating another  inversion  of  the  system  of  equations.  One  such  simulation  has  been  made 
with  the  HIM  code,  and  the  cost  was  increased  by  a factor  of  2.5 

Such  corrective  action  merely  solves  the  problem  of  large  rotation,  it  is  not  sufficient 
when  the  strains  themselves  are  large.  In  this  event  the  second-order  terms  must  be  included 
to  avoid  inaccuracies.  The  shear  strain  values  given  in  I aides  7-1  through  74  are  large,  indi- 
cating the  associated  linear  analyses  to  be  more  accurate  for  simulated  loads  scaled  down 
from  1 .000  pounds.  I o solve  the  problem  of  large  strain  in  the  brain  a nonlinear,  nearly 
incompressible,  finite  element  formulation  is  required  I his  element  has  bean  developed 
by  Hughes,  ct  al  ' * Some  numerical  experimentation  with  the  two-dimensional  version  of 
this  element  is  presented  in  the  following. 

\ single  quadrilateral  element  was  pinned  at  one  corner  and  assigned  an  initial  velocitv 
at  the  other  three  cornets  in  such  a way  as  to  induce  large  rotational  excursion  about  the 
pinned  corner  Because  of  the  centrifugal  forces  setup,  internal  strains  should  develop  in 
the  element  However,  with  the  linear  version  of  the  quadrilateral  element,  tin  computed 
strain  in  tlu  element  is  everywhere  /<•  ro,  and  the  attending  erroneous  dilation,  as  described 
above,  occurs  as  shown  in  I igurc  8 7 \V ith  the  nonlinear  element  the  results  are  quite  dil 
fereiit  tor  the  same  problem.  As  shown  in  f igure  8 8,  the  computed  l.agrangian  strain  I is 
non/cm  and  the  distortion  ol  the  square  element  is  consistent  with  the  centrifugal  force 
No  artificial  dilation  occurs  with  large  rotations  in  this  instance. 
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9 CONCLUSION'S 


A fully  three-dimensional  head  injury  model  computer  code,  referred  to  as  the  HIM 
code,  was  presented.  Its  development  from  preliminary  one-  and  two-dimensional  models 
through  to  the  present  completed  model  was  predicated  on  the  importance  of  recognizable 
geom  try  in  head  injury  modeling.  Time  histories  of  displacement  stress  and  strain  are  com- 
puted throughout  the  skull-brain  continuum  for  arbitrary  direct  impact  loads. 

The  model  is  applicable  to  elosed  brain  injury  exclusive  of  neck -related  trauma  as  well 
as  skull  bone  injury.  The  advantage  of  the  finite  element  approach  is  that  is  provides  a gen- 
eralized capability  for  studying  head  injury.  The  HIM  code  has  been  designed  with  suffi- 
cient flexibility  to  accommodate  new  information  as  it  becomes  available  regarding  geom- 
etry and  material  properties.  A wide  range  of  human  and  primate  skulls  can  be  easily  and 
automatically  discretized  to  assist  in  almost  any  head-impact  study.  The  accuracy  achieved 
in  such  studies  is  significantly  enhanced  by  employing  recognizable  geometry.  Boundary 
conditions  and  impact  loads  require  only  a minimum  of  idealization.  Any  or  all  parts  of  the 
skull-brain  system  can  be  idealized  as  being  either  viscoelastic  or  linear  elastic. 

A direct  integration  technique  is  used  to  integrate  the  equations  of  motion  in  the  HIM 
code.  I hc  primary  advantage  of  the  method  is  that  it  can  integrate  those  vibration  modes 
which  are  most  prominent  in  a response  due  to  an  arbitrary  impact  force,  whether  it  be  of 
short  or  long  duration.  Control  is  provided  by  specification  of  a time  step  size  consistent 
with  the  frequency  content  of  loading.  Besides  this  flexibility,  the  direct  integration 
method  can  be  readily  employed  in  an  extension  to  nonlinear  analyses. 

Results  from  a two-dimensional  axisvmmetric  finite  element  model  compared  favorably 
with  results  from  a similar  finite  difference  model.  In  those  instances  when  an  axisvmmetric 
fluid-filled  shell  may  be  considered  a useful  head  injury  model,  the  finite  element  method  is 
shown  to  be  competitive  with  other  numerical  and  analytical  approaches.  When  dealing 
w ith  load  durations  exceeding  1 ms,  it  is  perhaps  more  attractive  than  the  finite  difference 
method  and  in  most  cases  always  offers  more  flexibility  than  analytical  approaches. 

A two-dimensional  plane  strain  model  demonstrated  the  importance  of  simulating  more 
complex  geometry.  Skull  bending  and  intracranial  pressure  distributions  were  predicted 
more  realistically  than  in  previous  modeling  efforts.  Potential  sites  of  linear  skull  fracture 
and  closed  brain  injure  were  indicated  relative  to  anatomical  features.  Furthernvote,  the 
sensitivity  of  these  locations  to  the  characteristics  of  loading  and  boundary  conditions  w as 
clearly  evident 

A modification  of  the  isoparametric  element  formulations  for  simulating  brain  matter 
produced  a response  representative  of  clinically  observed  intracranial  response.  I he  clement 
was  obtained  by  a reduced  form  of  numerical  integration  termed  one-point  Gaussian  quad- 
rature. This  procedure  is  necessary  to  effect  a conventional  finite  element  analysis  of  nearlv 
incompressible  materials.  In  this  way  the  brain  is  idealized  as  a nearly  incompressible  mater- 
ial consistent  with  available  experimental  data  on  brain  material  properties. 

Loup  and  contrccoup  intracranial  response  are  very  evident  in  simulations  w ith  the 
model.  Inclusion  of  the  foramen  magnum  causes  a posterior  shilt  in  the  null  point  of  the 
otherwise  symmetrical,  anterior  posterior,  pressure  distribution.  This  behavior  is  consistent 
with  previous  experimental  observations,  f urther,  it  is  an  explanation  of  the  asymmetry  in 
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statistical  distribution,  noted  in  the  literature,  of  contrccoup  and  coup  injury  occurring  as 
a result  of  frontal  and  occipital  impacts. 

1‘hc  brain  is  simulated  as  being  isolated  from  the  cranial  wall  by  including  a compressible 
representation  of  the  subarachnoid  space.  In  this  way  intracranial  volume  shifts  and  the 
brain  slosh  mechanisms  that  have  been  clinically  observed  can  also  be  observed  in  the 
model’s  response. 

Automated  mesh  generators  are  the  only  practical  means  for  developing  discretizations 
tor  head  injur)'  models  which  depict  recognizable  geometry.  They  cannot  avoid  all  the 
laborious  effort  of  extracting  measured  data  because  a minimum  of  this  information  is  nec- 
essary. Nevertheless,  they  are  to  be  preferred  over  completely  manual  efforts  at  generating 
the  discretizations.  Automatic  mesh  generators  result  in  skull-brain  discretizations  which 
are  uniform  as  well  as  consistent  with  the  specified  geometry. 

Techniques  for  efficient  equation  solving,  including  optimal  renumbering  of  the  system 
of  motion  equations,  arc  necessary  for  an  economically  flexible  head-injury  model.  Two 
basic  discretizations  are  provided  the  user  in  the  HIM  code:  a fine  and  a coarse  mesh.  F.ach 
allows  for  taking  full  advantage  of  symmetrical  conditions  when  they  exist.  The  present 
linear  HIM  code  costs  approximately  $200  per  simulation,  including  apprximatcly  $30  of 
data  post-processing.  These  costs  are  based  upon  present  commercial  computer  rates  and 
involve  anywhere  from  10  to  15  time  steps  with  the  fine  discretization. 

The  HIM  code  has  undergone  extensive  numerical  testing,  and  all  known  bugs  have  been 
| corrected.  With  the  exception,  perhaps,  of  simulating  the  brain  material,  the  HIM  code  has 

been  constructed  with  standard  dynamic,  finite  element  techniques.  A qualitative  validation 
has  been  established  through  comparisons  w'ith  theory  for  simple  models,  through  compari- 
son with  solutions  from  methods  other  than  finite  element,  and  through  observing  the  rea- 
sonableness of  both  the  skull  bone  response  and  the  coup  and  contrccoup  response  in  the 
brain  Attempts  at  quantitative  validation  have  not  proved  successful.  I he  present  skull 
bone  finite  element  characterization  is  found  to  be  too  stiff.  To  correct  the  problem  an 
eight-node  incompatible  element  has  been  developed  which  will  allow  more  bending  but, 
as  yet,  has  not  been  incorporated  and  tried.  Comparison  of  measured  and  computed  intra- 
cranial pressures  are  difficult  but  were  attempted  anyway  in  this  study.  Measured  epidural 
pressures  in  anesthetized  rhesus  monkeys  were  on  the  order  of  10  psi  while  computed  intra- 
cranial pressures  of  primate  simulations  were  on  the  order  of  37  psi  While  this  is  only  with- 
in an  order  of  magnitude  agreement,  previous  model  attempts  were  known  to  have  com- 
puted pressure  in  excess  of  100  and  even  1,000  psi.  Future  improvements  in  both  experi- 
mental and  modeling  techniques  may  well  establish  a satisfactory  quantitative  validation. 

Simulations  with  the  HIM  code  have  assisted  in  defining  the  li/nits  of  its  applicability 
i imposed  by  the  assumptions  of  linearity.  Rigid  body  rotations  in  excess  of  10  degrees  may 

require  nonlinear  strain  terms  if  the  response  of  interest  has  not  already  subsided.  Shear 
strains  were  predicted  which  are  in  the  realm  of  large  strain  analysis  anti  indicate  a require 
t ment  for  nonlinear  formulations  independent  of  large  rotations.  Ffforts  are  presently 

underway  to  effect  a nonlinear  analysis.  Meanwhile  the  linear  HIM  code  predictions  can 
yield  useful  insight  into  the  mechanical  causes  of  skull  and  brain  injury. 

Qualitative  information  obtained  from  the  linear  HIM  code  can  provide  useful  insight 
into  the  mechanical-clinical  relationships  of  head  injury.  Though  this  report,  in  the  main, 
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describes  the  development  of  a head  injury  analysis  tool,  some  results  derived  during  the 
course  of  the  project  do  fall  under  this  category  and  are  discussed  in  the  following  para- 
graphs. 

1 he  skull  bone  stress  response  induced  by  an  occipital  impact  suggests  that  impending 
skull  fracture  would  occur  along  a path  between  the  point  of  loading  and  the  foramen  mag- 
num. While  the  skull  is  thick  in  this  region,  the  foramen  magnum  presents  an  area  of  stress 
concentration. 

I lie  relative  displacement  between  the  skull  and  brain  observed  in  the  predicted  re- 
sponse has  also  been  observed  in  tests  involving  primates  with  implanted  lucite  calvariums. 
Th is  motion  appears  beneficial  in  alleviating  somewhat  the  stress  in  the  brain.  But  the  bene- 
ficial effect  is  true  only  up  to  a point.  Strain  develops  at  the  brain’s  surface  when  further 
relative  motion  is  prevented  by  either  geometrical  or  mechanical  constraints  at  the  skull- 
brain  interface.  The  tendency  for  brain  elements  nearest  this  interface  to  deform  more  than 
those  in  the  midbrain  area  is  very  apparent  in  the  graphics  of  the  deformed  mesh. 

The  boundary  condition  specification  has  been  shown  to  influence  significantly  the 
computed  intracranial  response  by  a comparison  of  a restrained  skull,  a partially’  restrained 
skull,  and  an  unrestrained  skull,  all  subjected  to  the  same  impact  load.  Predicted  intra- 
cranial stress  distributions  develop  null  pressure  regions  which  tend  to  locate  above  the 
foramen  magnum  independent  of  the  direction  of  impact.  Thus,  for  unrestrained  impacts, 
significant  contrecoup  tensile  stresses  will  exist  for  impacts  to  the  occiput  and  significant 
| coup  compression  stresses  will  exist  for  impacts  to  the  frontal  bone.  The  contours  at  the 

anterior  and  posterior  brain  surfaces  locate  closely  the  sites  of  potential  coup  and  contre- 
coup damage  in  the  midsagittal  plane. 

Shear  strain  predictions  for  direct  impact  loads  have  shown  the  brainstem  region  to  be 
most  vulnerable.  Besides  the  brainstem  region,  the  outer  surfaces  of  the  brain  developed 
shear  strains.  1 he  midbrain  region  evidenced  little  or  no  shear  strain.  This  response  is  true 
also  of  the  indirect  impact  load  case  except  that  the  maximum  shear  strains  occurred  at  the 
brain  surface  as  opposed  to  the  brainstem  region. 

\ The  results  are  generally  independent  of  the  restraint  conditions  at  the  base  of  the  skull 

and  thus  the  neck  is  not  expected  to  modify  them  greatly. 
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10.  RECOMMENDATIONS 


Recommendations  can  be  divided  into  two  categories.  The  first  category  pertains  to  the 
linear  I IIM  code  and  the  second  category  pertains  to  the  development  of  a nonlinear  head 
injure  model. 

f irst,  with  regard  to  improving  the  linear  HIM  code  and  its  validation  with  experimental 
data,  three  specific  recommendations  are  made. 

1 . The  linear  isoparametric  element,  eight-node  brick  used  to  model  the  skull  bone 
structure  should  be  replaced  with  an  incompatible  finite  element  formulation  to  allow  more 
effective  modeling  of  the  bending  stress  gradient  across  the  simulated  three-layered  cranial 
thickness.  The  skull  clement  should  remain  three-dimensional,  as  is  now  the  case,  as  op- 
posed to  a two-dimensional  shell  element. 

2.  Employ  reduced  integration  (one-point  quadrature)  wifh  the  dilatational  terms  of 
the  brain  element  stiffness  matrix  only,  while  using  the  normal  quadrature  rule  for  the  dis- 
tortional  energy’  terms.  This  should  preserve  the  present,  nearly  incompressible,  formula- 
tions for  the  brain  element  while  reducing  the  strain  magnitudes,  which  are  typically  very’ 
large  in  the  present  HIM  code  computations. 

3.  Improvement  of  the  HIM  code  running  efficiency  can  be  made  by  the  following: 

(a)  Incorporating  the  bandwidth  minimizing  routines  directly  into  the  HIM  code 

so  that  their  use  is  transparent  to  the  user. 

I (b)  A more  efficient  procedure  can  be  adapted  for  solving  the  equilibrium  equation 

if  special  hardware  features  available  on  the  CDC  6600  are  employed. 

(c)  Direct  access  I/O  should  be  used  in  processing  the  right-hand  side  vectors.  I his 
would  be  especially  useful  for  dynamic  analyses  of  multiload  vectors. 

(d)  The  inner  product  calculations  involved  in  backward  reduction  could  be  done 
more  efficiently.  The  same  techniques  noted  in  (b)  above  should  be  adapted. 

■ (e)  I/O  methods  currently  employed  in  the  program  could  be  reworked  advan- 

tageously. 

(0  The  program  should  be  overlayed. 

With  these  improvements  it  is  estimated  that  the  running  costs  could  be  reduced  by 
50%.  Further  recommendations  would  include  the  following. 

1 Injury  or  “failure”  criteria  for  the  simulated  structures  of  the  brain  and  skull  should 
be  developed  and  incorporated  into  the  HIM  code;  then  the  computed  results  summarized 
into  a printed  description  of  injury  hazard  for  each  HIM  code  simulation. 

! 2.  Additional  parameter  studies  should  be  conducted  for  the  purpose  of  generating  a 

head  injury  data  base  which  would  facilitate  directly  the  preparation  of  safety  standards. 

3.  Preprocessor  and  postprocessor  integration  into  the  HIM  code  would  greatly  en- 
hance use  of  the  code  and  analysis  of  the  data.  These  routines  have  been  completed,  but 
remain  separate  from  the  HIM  code. 

Recommendations  pertaining  to  the  nonlinear  head  injury  model  are  less  specific  than 
the  above  recommendations.  A nonlinear  computer  code  similar  to  the  basic  HIM  code 
exists,  but  the  HIM  code  input  module  (geometry,  etc.)  must  be  incorporated  and  checked 
out  first  before  more  detailed  and  specific  recommendations  can  be  made  regarding  a non- 

t linear  code.  Assuming  this  is  accomplished,  the  following  recommendations  can  be  made: 
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1 . Repeat  selected  linear  HIM  code  runs  with  the  nonlinear  model  to  discern  the  dif- 
ferences anil  further  define  the  limitations  of  the  less  expensive  linear  model. 

2.  With  the  large-rotation/large-strain  capability,  simulate  the  large  head  rotation  data 
typical  of  NAMR1.  test  sled  runs. 

3.  Investigate  appropriate  temporal  integration  techniques  to  eliminate  numerical 
noise  associated  with  the  use  of  the  Newmark  method  along  with  prescribed  displacement 
data  input. 

4 Develop  a capability  to  input  prescribed  acceleration  data  to  eliminate  any  error 
associated  with  numerical  smoothing  (double  integration)  of  measured  acceleration  data. 

5 Develop  a capability  to  synthesize  the  nonlinear  head  injury  model  code  and  the 
contact/impact  prediction  code.  In  this  way,  the  dynamic  impact  of  both  the  simulated 
skull  and  target  can  be  included  in  one  computer  run. 
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CLINICAL  DESCRIPTION  OF  HEAD  INJURY 


This  appendix  is  a reproduction  of  a report  submitted  in  March  1974  by  Richard  B. 
Small,  M.D.  of  the  University  of  Southern  California  Medical  Center  in  fulfillment  of  CEL 
Contract  No.  N62  399-7  3-C 4)006. 
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INTRODUCTION 


Historically  man  has  long  recognized  the  serious  nature  of  head  injuries.  F arly  human 
societies  developed  mechanisms  for  dealing  with  head  injury  as  evidenced  by  the  develop- 
ment of  helmets  for  protection  of  the  head  in  warfare.  Some  of  these  societies  also  prac- 
ticed trepanation,  however  its  use  cannot  be  related  to  head  injury.  Later  cultures  have 
recorded  rather  sophisticated  observations  regarding  head  injury.  In  the  Edwin  Smith  sur- 
gical papyrus  (circa  1600  BC)  several  cases  of  head  injury  arc  presented  along  with  the  ob- 
servation of  contralateral  weakness  anil  the  poor  prognosis  associated  with  feeble  pulse  and 
fever.  In  the  Hippocratic  writings  there  arc  a number  of  observations  referring  to  head  in- 
jury and  the  first  treatise  devoted  solely  to  the  subject.  Included  are  recommendations  for 
elevation  of  depressed  fractures  and  trepanation.  Attention  was  not  directed  to  the  neuro- 
logical condition  of  the  patient  until  much  later.  Delineation  of  the  presently  recognized 
syndromes  comprising  the  post -traumatic  encephalopathies  was  not  accomplished  until  the 
late  18th  Century'. 

Commotio  Cerebri.  The  symptomatology  of  concussion  of  the  brain  was  presented  by 
Lanfranci  of  Milan  in  his  Chirgia  Magna  (1296  AD).  This  was  presented  as  a classic  rather 
than  an  original  description.  Essentially  the  same  description  can  be  found  in  the  writings 
of  Ambrose  Pare  (1561),  a number  of  18th  Century  writers  including  Littre  (1705),  Jean- 
l ouis  Pettit  (1774)  and  Benjamin  Bell  (1783),  anil  recent  reviews  including  Courvillc  (1951 ) 
and  Ward  ( 1966). 

Contusio  Cerebri.  Concerning  cerebral  contusions  Pare  (1561)  recognized  that  in  some- 
injuries  “by  the  violence  of  the  blow  the  veins  and  arteries  may  be  broken.”  Shortly  there- 
after l allopius  ( 1569)  described  similar  findings  contralateral  to  the  impact.  The  signifi- 
cance of  the  observations  was  well  recognized  prior  to  1766  when  the  Academic  Royale  dc 
Paris  offered  a prize  for  “ctablir  la  theorie  dea  contracoups  dans  les  lesions  de  la  tete  et  les 
consequences  qu’on  peut  cn  titer.”  T his  subject  has  continued  to  be  of  considerable 
interest. 

Compressio  Cerebri.  It  is  possible  (Courvilie,  1951 ) that  Celsus  recognized  that  intra- 
cranial hemorrhage  or  sanguanous  effusion  resulted  from  head  injury.  Pare  described  sub- 
dural hematoma  and  suggested  tearing  of  veins  as  its  etiology  (Flamm,  1972).  The  first 
definite  distinction  of  the  symptomatology'  of  compression  from  concussion,  however,  is 
found  in  l.e  Dran  (1731). 

Nevertheless  considerable  confusion  has  persisted  to  the  present  time  concerning  the 
interrelationship  of  cerebral  concussion,  contusion  and  compression.  T his  is  manifested  in 
the  concept,  which  developed  in  the  middle  of  the  19th  Century  and  persists  today,  of  these 
entities  as  a continuum  particularly  in  clinical  evaluation.  Anatomically  the  lesions  are  dis- 
tinctly different.  Biomechanically  their  linear  relationship  is  not  substantiated  by  present 
evidence. 

Experimental  studies  of  head  injuries  were  carried  out  during  the  19th  Century  by  (lama 
(1835)  and  others,  but  well  controlled  studies  were  not  performed  until  earlier  in  the  pres- 
ent century.  Basic  to  modern  studies  was  the  ncuropathological  work  of  Ramon  y Cajal. 

T his  work  was  extended  to  human  material  by  Spatz  and  co-workers  in  Germany  and  Rand 
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and  Courv’lle  in  this  country.  Clinical  and  experimental  studies  by  Denny-Brown  and  Bus- 
sell (1932,  1940,  and  1945)  stimulated  the  investigation  of  the  clinico-pathologico- 
tnechanical  correlation  of  head  injury,  which  is  the  ultimate  objective  of  this  present 
project. 

Theoretical  Considerations.  The  mechanism  by  which  the  impact  of  a closed  head  in- 
jury produces  damage  and  dysfunction  in  the  central  nervous  system  has  been  a subject  of 
interest  to  clinicians  and  pathologists  for  some  time.  That  there  is  no  concensus  of  opinion 
is  reflected  in  a rather  extensive  literature.  Early  vascular  and  humoral  theories  have  been 
discarded,  and  there  is  general  acceptance  that  mechanical  forces  produce  intracranial  in- 
juries. The  oldest  mechanical  theory  is  that  of  vibration  which  was  the  consensus  of  the 
Parisian  Royal  Academy  of  Surgery.  Another  supporter  was  Gama  who  developed  the  first 
experimental  model  of  head  injury  in  1835.  Me  imbedded  black  threads  in  a gellatin  filled 
flask.  To-and-fro  movement  of  the  threads  was  noted  when  the  flask  was  struck.  The  theor- 
ies derived  from  this  model  included  vibration  of  the  skull,  pressure  gradients  in  the  brain 
and  propagation  of  waves  of  force  or  compression  through  the  brain.  Experimental  mea- 
surements have  tended  to  discredit  these  theories  as  the  pressure  changes  measured  were 
propagated  too  slowly  to  explain  effects  observed  with  intermediate  high  speed  cameras 
(1000  - 4000  frames  per  second).  Such  studies  have  been  performed  by  Pudcn/.  and 
Sheldon  ( 1944),  Ommaya  ( 1969)  and  Gosch  et  al  ( 1970).  They  have  confirmed  the  move- 
ment of  the  brain  relative  to  the  skull  postulated  by  Russell  (1932)  and  Molbourn  (1943- 
44).  Molbourn  is  the  first  of  a number  of  physicists  who  have  contributed  a great  deal  to 
the  study  of  the  biomechanics  of  head  injury.  Mis  theory  is  based  on  shear  forces  developed 
within  the  cranium  by  rotational  acceleration  of  the  head.  Such  forces  arc  also  invoked  in 
the  explanation  of  the  axonal  injuries  described  by  Strich  and  others.  Many  authors  have 
suggested  less  severe  and  reversible  damage  by  shear  forces  to  fiber  systems  in  mild  and  mod- 
erate degrees  of  post  traumatic  encephalopathy.  Rotational  acceleration  has  also  been  used 
with  some  success  in  describing  the  distribution  of  contusions,  further  evidence  of  the 
existence  of  shear  forces  in  head  injury  is  derived  from  high  speed  camera  studies  of  photo- 
elastic models. 

Another  physical  scientist  suggested  a theory  based  on  a simple  model  not  unlike 
Gama’s.  Gross’  model  was  a fluid  filled  flask  in  which  intermediate  high  speed  photography 
revealed  cavitation  following  impact.  Collapse  of  these  cavities  is  postulated  to  produce 
high  local  accelerations  leading  to  shear  forces,  further  development  of  this  theory  has 
been  carried  out  by  Unterharnscheidt  and  Sellicr  (1966)  and  by  Kopccky  and  Rippergcr 
(1969).  In  their  model  the  negative  pressures  produced  by  acceleration  and  elastic  defor- 
mation of  the  skull  do  not  necessarily  drop  so  low  as  to  produce  cavitation,  however,  they 
produce  pressure  differentials  across  cell  membranes  sufficient  to  cause  disruption  of  their 
integr’ty. 

A divergent  mechanism  of  injury  is  invoked  very  successfully  by  Lindenberg  (1960)  to 
explain  a selected  group  of  fatal  human  head  injuries.  Thus  elastic  deformation  of  the  skull 
and  translational  acceleration  or  “mass  movement”  of  the  brain  in  the  “line  of  impact”  is 
shown  to  be  reasonable  mechanism. 

The  brain  is  not  rigidly  supported  within  its  own  confines,  but  is  tethered  somewhat  at 
its  base  by  the  exit  of  cranial  nerves  and  the  spinal  cord  from  the  cranium.  Bridging  veins 
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from  the  dural  venous  sinuses  to  the  cortex  add  another  minor  stabilizing  force.  The  falx 
and  the  tentorium  lend  more  major  support.  Yet  the  brain  is  a semisolid  medium  with 
surrounding  cerebrospinal  fluid.  Blows  to  the  head  cause  significant  shifts  of  brain  tissue 
within  the  cranium.  Let  us  consider  the  intracranial  compartment. 

The  skull  is  a rigid  cavity  with  a capacity  of  1400  - 1500  cc’s,  an  irregular  base,  and  a 
smooth  vault  which  is  composed  of  an  inner  and  an  outer  table  of  bone  with  innerposed 
cancellous  bone.  This  arrangement  provides  greater  strength  with  an  economy  of  weight, 
as  well  as  a certain  degree  of  insulation.  The  cranial  vault  is  not  uniform  in  its  thickness, 
but  varies  from  2 mm.  in  the  temporal  region  to  greater  than  1 cm.  in  the  mid-frontal  and 
parietal  occipital  regions.  Irregularities  of  the  base,  particularly  the  sphenoid  ridges  and  the 
orbital  plate  areas,  greatly  modify  the  formation  of  fractures,  and  the  overall  planes  of 
energy  transfer  through  the  skull  and  thus  through  the  brain. 

The  base  of  the  skull  is  a membranous  structure  with  a cancellous  or  diploic  layer  only 
in  certain  areas.  The  three  main  divisions  of  the  base  are  the  anterior,  middle,  and  the  pos- 
terior fossae.  The  anterior  fossa  is  formed  by  the  horizontal  or  orbital  portion  of  the  frontal 
bone  laterally,  the  lesser  wing  of  the  sphenoid  and  the  body  of  the  sphenoid  posteriorly,  and 
the  cribiform  plate  of  the  ethmoid  in  the  midlinc.  Projecting  upward  from  the  midline  of 
the  cribiform  plate  in  its  forward  portion  is  the  crista  galli,  to  which  is  attached  the  falx 
cerebri.  The  middle  fossa  is  formed  by  the  greater  wing  of  the  sphenoid  (containing  the 
pituitary  fossa)  in  the  middle,  and  squama  and  the  anterior  aspect  of  the  petrous  portion  of 
the  temporal  bone  laterally  and  posteriorly.  The  anterior  and  middle  fossae  arc  separated 
by  the  sphenoid  ridge,  with  its  sharp,  knife-like  edge  fitting  against  the  junction  of  the 
frontal  and  temporal  poles  of  the  cerebral  hemispheres.  The  posterior  fossa  is  formed  by 
the  posterior  aspect  of  the  petrous  portion  of  the  temporal  bone  anteriorly  and  the  occipital 
bone  posteriorly.  The  foramen  magnum  is  in  the  center  of  the  posterior  fossa,  with  the 
basilar  process  of  the  occipital  bone  and  the  posterior  part  of  the  body  of  the  sphenoid 
forming  the  base  of  the  skull  in  front  of  the  foramen  magnum.  To  further  complicate  the 
area,  numerous  perforations  which  provide  avenues  of  egress  for  cranial  nerves  arc  distri- 
buted acrosss  the  base.  The  main  vascular  supply  to  the  intracranial  contents  also  pene- 
trates the  skull  via  this  area.  Fractures  occurring  through  the  base  may  produce  injury  to 
any  of  the  cranial  nerves  or  vascular  structures  which  arc  intimately  associated  with  the 
area.  These  injuries  are  produced  cither  by  direct  shearing  effects,  or  by  the  accumulation 
of  hemorrhage  within  a foramina  with  subsequent  compromise  of  space  and  compression. 

In  addition,  fractures  through  this  area  may  produce  dural  tear  with  cerebrospinal  fluid  leak 
via  a route  of  cranial  nerve  egress. 

Linear  fracture  lines  tend  to  occur  radially  from  a center  of  inbending  or  deformation 
secondary  to  an  energy  transfer  to  the  skull.  These  fractures  are  modified  by  areas  of  bony 
strength.  Fracture  lines  are  predetermined  by  the  site  of  impact  as  well  as  the  strength  of 
the  skull.  Most  often: 

1 . Blows  frontally  cause  linear  fractures  into  the  orbital  roofs; 

2.  Vertical  blows  cause  radiating  fractures; 

3.  Occipital  blows  cause  fractures  running  toward  the  base  and  foramen  magnum; 

4.  Parietal  blows  cause  fractures  directed  toward  the  temporal  region. 
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Depressed  fractures  are  usually  caused  by  a relatively  sharp  point  of  contact  with  lacer- 
ation of  scalp  being  common.  With  massive  quantities  of  energy  transfer,  large  flat  objects 
may  cause  extensive  depression  of  bone. 

In  considering  the  intracranial  volume,  dural  and  meningeal  supporting  structures  may 
be  considered  as  inconsequential.  Appropriate  values  for  the  remaining  intracranial  contents 
include: 

1.  Brain  88.8%  (65%  H20); 

2.  Blood  2,3%; 

3.  Cerebrospinal  fluid  8.9%. 

The  skull  is  rigid  and  allows  only  minimal  expansion  through  the  foramen  magnum  and 
the  smaller  foramina  that  transmit  blood  vessels  and  nerves.  In  addition,  dural  reflections 
divide  the  intracranial  cavity  into  fossae  that  normally  protect  the  brain  against  excessive 
movement  but  limit  the  amount  of  compensary  shift  and  displacement  that  can  develop  in 
response  to  abnormal  conditions. 

In  this  regard,  the  tentorium  cerebelli  divides  the  cranial  vault  into  anterior  and  poster- 
ior fossae.  The  rigid  dural  extension  leads  posteriorly  from  the  petrous  ridges  and  anterior 
clinoid  processes,  sloping  downward  and  laterally  from  its  medial  edge  to  attach  to  the  occi- 
pital bone  along  the  line  of  the  lateral  sinus.  Extending  posteriorly  into  the  center  of  the 
tentorium  from  the  posterior  clinoid  processes  is  a large  semioval  opening,  the  incisura  or 
tentorial  notch,  whose  diameters  are  usually  between  50  and  70  mm.  in  the  fronto-occipital 
axis  and  25  ami  40  mm.  in  the  interparietal  axis. 

The  temporal  lobes  rest  on  the  tentorial  incisura  and  their  medial  surfaces  protrude  over 
its  edge,  so  that  3 to  4 mm.  of  the  medial  anterior  portion  of  the  temporal  unous  bulges 
into  the  notch.  A small  segment  of  the  hippocampal  gyrus  also  overhangs  the  edge  and  be- 
comes more  narrow  posteriorly. 

It  is  important  to  retain  the  fact  that  changes  in  the  relationships  between  the  tentorial 
incisura  and  its  surroundings  explain  most  of  the  complications  in  supratentorial  volume 
shifts.  The  midbrain  occupies  the  anterior  portion  of  the  notch.  The  cerebellum  lies  in 
juxtaposition  to  the  dorsum  of  the  midbrain  and  fills  the  posterior  portion  of  the  notch. 
Ventral  to  the  brainstem  lies  the  basilar  artery,  which  divides  into  two  diverging  posterior 
cerebral  arteries  just  caudal  to  the  incisura.  I'.ach  posterior  cerebral  artery  crosses  the  ocu- 
lomotor nerve  that  emerges  caudal  to  it.  The  artery  then  circles  the  homolateral  cerebral 
peduncle  and  the  adjacent  lateral  midbrain  and  reaches  the  vcntal  surface  of  the  hippo- 
campal gyrus  of  the  temporal  lobe,  where  it  crosses  the  tentorial  edge  and  proceeds  toward 
the  occipital  lobe. 

The  oculomotor  nerves  penetrate  the  medial-basal  surface  of  each  cerebral  peduncle 
caudal  to  the  tentorium.  They  traverse  the  basal  cistern  (tentorial  gap),  initially  passing  be- 
tween the  superior  cerebellar  and  posterior  cerebral  arteries,  and  then  passing  in  close  prox- 
imity to  each  temporal  lobe  uncus  at  the  point  where  the  uncus  overhangs  the  lateral  in 
cisural  edge,  and  finally  passing  over  the  petroclinoid  ligaments  to  enter  the  cavernous  sinus. 

(ionsidering  that  the  brain  may  move  within  the  skull,  it  is  evident  that  two  types  of 
injury  may  be  produced  by  a blow  to  the  skull: 

1 . If  energy  is  sufficient  at  the  point  of  impact  the  brain  will  be  damaged  to  a degree 
compatible  with  the  force.  This  is  called  a COUP  type  injury. 
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2.  The  delivered  energy  will  cause  the  brain  to  shift  in  a direction  dictated  by  the  force- 
vector  of  injury,  causing  it  to  rebound  off  the  inner  cranial  wall,  and  inflicting  injury  pro- 
portional to  the  force  and  related  to  the  configuration  of  the  particular  section  of  the  cra- 
nial wall,  a CONTRECOUP  type  injury. 

The  importance  of  intracranial  volume  shifts  in  the  production  of  brain  injury  cannot  be 
overestimated.  Vascular  and  neural  injury  secondary  to  these  displacements  are  often  re- 
sponsible for  injuries  which  produce  irreconcilable  results.  Mesial  shift  of  the  lateral  hemi- 
sphere may  force  the  cingulate  gyrus  under  the  falx  cerebri.  Downward  displacement  of  the 
hemispheres  and  basal  nuclei  may  compress  and  displace  the  dienccphalon  and  the  adjoining 
midbrain  rostrocaudally  through  the  tentorial  incisura.  In  such  transgressions  the  pituitary 
stalk  may  be  partially  avulsed  and  marked  changes  may  occur  in  the  diencephalon  and  brain- 
stem. More  posterior  mesial  shift  may  result  in  compression  of  the  midbrain  against  either 
the  free  edge  of  the  tentorium  of  some  other  resistant  structure  such  as  the  petroclinoid 
ligament.  Herniation  at  the  incisura  also  has  direct  effects  upon  the  brain  stem  and  upon  the 
delicate  vascular  penetrations  of  the  structures  as  well.  Stretching  of  the  perforating 
branches  of  the  basilar  artery  occurs  because  the  artery  is  tethered  to  the  circle  of  Willis  and 
cannot  shift  caudally.  Thus,  infarction  and  hemorrhage  are  produced  within  the  brain  stem 
to  further  complicate  existing  injury. 
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figure  A-l.  The  Cranium  liemisection  of  the  skull. 
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I’lgurc  A 7 I he  Meninges  and  Blood  Vessels  the  vessels  over  the  medial  surface. 
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figure  A 9 I Ik-  Meninges  and  Blood  Vessels  the  posterior  fossa  and  cervical  cord. 
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